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SECTION  1 


OBJECTIVES 


An  analytical-numerical  study  was  pursued  by  the  present  investigators  under 
AFOSR  Grant  No.  F  49620-92-J-0292  between  May  1992  and  December  1995.  The 
primary  objectives  were:  (i)  to  characterize  the  unsteady  separation  in  terms  of 
vortex  dynamics  that  leads  to  dynamic  stall  of  an  airfoil  and  a  rectangular  wing  in 
constant  pitch-rate  motion;  (ii)  to  identify  the  physical  mechanism  responsible  for 
moment  stall,  and  formulate  an  active-control  strategy  to  manage  stall.  For  this 
program,  one  additional  graduate  student's  effort  was  made  possible  through  the 
AASERT  Grant  No.  F  49620-93-1-0398,  between  June  1,  1993  and  May  31,  1996. 
This  grant  was  extended,  at  no  cost,  for  one  additional  year,  to  May  3 1,  1997.  For 
this  AASERT  grant,  the  major  objective  was  to  improve  the  analysis  tools  as  well  as 
assure  the  level  of  accuracy  and  efficiency  that  can  be  realized  in  the  study  of  physics 
of  the  dynamic-stall  and  related  unsteady  phenomena  through  the  development  of  the 
LES/DNS  methodology. 

The  effort  on  the  AASERT  project  was  started  with  Ms.  Susan  Beltz  who 
studied  unsteady  flows  using  a  vorticity-based  formulation.  She  examined,  in  detail, 
the  surface  and  outflow  boundary  condition  for  the  unsteady  formulations,  and  wrote 
a  Master's  Thesis  in  August  1994.  Due  to  an  exciting  offer  from  a  high-tech  software 
company,  run  by  my  previous  students  who  had  worked  on  the  AFOSR  project, 
Susan  Beltz  decided  to  temporarily  bifurcate  from  pursuing  a  doctoral  degree 
program.  Keith  Blodgett  a  doctoral  student,  was  appointed  on  this  project  in  place 
of  Susan  Beltz.  Keith  worked  on  a  spectral  technique  in  generalized  coordinates,  also 
using  a  vorticity-based  formulation,  to  study  the  problem  of  receptivity  and  completed 
his  Doctoral  Dissertation  in  May  1995.  Subsequently,  Chris  Noll  and  Brad  Duncan, 
both  doctoral  students,  were  partially  supported  on  this  project.  They  will  complete 
their  doctoral  dissertations  by  March  1998.  Brad  Duncan  has  been  using  direct 
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numerical  simulation  methodology  and  a  spectral  technique,  in  the  area  of  early  stages 
of  transition,  before  the  nonlinear  breakdown  occurs.  On  the  other  hand,  Chris  Noll 
has  been  studying  the  effect  of  compressibility  on  the  flow  past  maneuvering  bodies. 
The  last  two  studies  are  heavily  involved  in  the  new  paradigm  of  Object-Oriented 
Programming  (OOP).  Chris  Noll  has  been  improving  the  OOP,  and  has  created  an 
Object  Oriented  Numerics  (OON)  that  not  only  facilitates  large-scale  programming 
requiring  complex  software,  but  also  has  been  able  to  retain  high  code  efficiency. 

Through  GPA,  both  Chris  Noll  and  Brad  Duncan  have  been  able  to  develop 
software  that  uses  all  the  available  tools  of  high-performance  computing,  and  have 
also  developed  parallel  versions  of  this  software  for  more  than  one  high-performance 
platform. 


Each  of  these  four  studies  are  briefly  described  next. 
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SECTION  2 


DESCRIPTION  OF  SIGNIFICANT  ACCOMPLISHMENTS 


All  of  the  areas  of  research  pursued  and  the  progress,  as  well  as  specific 
achievements,  made  in  these  studies  during  the  four-year  grant  period  are  briefly 
summarized  in  this  section. 

Since  the  technical  discussion  pertaining  to  progress  made  for  some  individual 
objectives  is  quite  long,  the  details  are  relegated  to  the  corresponding  reports  and 
papers  included  in  the  appendices  of  this  report. 

2A.  Improving  the  Outflow  Boundary  Condition  for  Vorticity-Based 
Formulation  -Susan  Beltz,  M.S.  Student 

A  vorticity-stream  function  circulation  ( (o,  )  formulation  was 

developed  based  upon  our  earlier  work,  which  focused  on  the  accurate 
implementation  of  the  far-field  boundary  condition  to  provide  a  unique  formulation. 
Although  the  study  considered  unsteady  incompressible  viscous  flows  past  Joukowski 
airfoils,  the  ability  of  the  formulation,  which  used  generalized  orthogonal  coordinates, 
to  incorporate  any  sequence  of  conformal  mappings  enables  a  range  of  flow 
configurations  to  be  studied.  The  uniqueness  in  the  formulation  has  been  realized 
through  the  incorporation  of  viscous  circulation  into  the  formulation.  The 
resulting  system  of  governing  equations  was  solved  numerically  using  the  Block 
Gaussian  Elimination  (BGE)  procedure  to  solve  the  integral  form  of  the  curl  equation, 
and  the  Alternating  Direction  Implicit  (ADI)  technique  of  Douglas  and  Gunn  to  solve 
the  vorticity-transport  equation,  with  an  integrative  procedure  used  to  ensure  pressure 
continuity  along  the  airfoil  surface.  Due  to  a  lack  of  available  data  pertaining  to 
Joukowski  airfoils,  the  validation  study  was  performed  using  the  circular  cylinder 
geometry,  for  which  there  exist  a  wealth  of  available  experimental  and  numerical  data. 
Specifically,  validation  was  performed  on  the  cylinder  geometry  using  early-time 
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results  for  two  different  values  of  the  Reynolds  Re  number.  Based  on  measurements 
of  the  length  and  width  of  the  symmetric  eddies  which  form  at  the  rear  of  the  cylinder, 
a  comparison  was  made  with  the  available  experimental  data.  Additionally,  the  grid 
independence  of  the  solutions  obtained  was  analyzed  using  the  Grid  Convergence 
Index  (GCI).  It  was  found  that,  for  the  case  of  Re=9500,  a  grid  size  of  (705x91)  is 
adequate,  based  on  the  GCI.  Thus,  the  present  analysis  was  validated  and  is  available 
for  application  to  a  class  of  unsteady  two-dimensional  incompressible  flows.  The 
results  of  this  study  are  discussed  in  detail  by  Beltz  and  K  Ghia  (1995). 


2B.  A  Spectral  Multidomain  Method  for  Studying  Receptivity 

-  Keith  Blodgett,Ph.D.  Student 

A  spectral  method  in  generalized  curvilinear  coordinates  was  developed  for 
investigating  the  early  stages  of  transition  to  turbulence  on  complex  geometries.  The 
earliest  stage  of  transition  is  receptivity,  the  process  by  which  small  disturbances, 
originating  in  the  free  stream  penetrate  the  boundary  layer  and  excite  the  boundary- 
layer  instabilities.  A  spectral  method  was  chosen  over  other  numerical  methods 
because  of  its  higher  accuracy  and  negligible  phase  errors,  both  of  which  are  deemed 
important  in  the  simulation  of  receptivity  problems.  The  spatial  directions  were 
discretized  with  a  spectral  collocation  method  using  the  Gauss-Lobatto  points,  with  a 
choice  of  either  Chebyshev  or  Legendre  polynomials.  A  semi-implicit  finite-difference 
scheme  was  utilized  for  the  temporal  discretization,  with  the  linear  terms  being  treated 
implicitly  and  the  nonlinear  terms  treated  explicitly.  The  resulting  algebraic  system  of 
equations  was  solved  iteratively  sing  a  preconditioned  truncated  conjugate  residual 
method.  Several  issues  were  addressed,  including  the  choice  of  formulation  of  the 
governing  equations,  the  accurate  implementation  of  the  no-slip  boundary  condition, 
the  choice  of  preconditioner,  and  the  solution  method.  A  multidomain  method  was 
also  incorporated  in  to  the  spectral  method  to  aid  in  the  resolution  of  the  disparate 
length  scales  of  the  problem.  In  effect,  it  decomposes  the  problem  into  a  set  of  sub¬ 
problems,  which  are  related  through  interface  conditions.  Two  interface  conditions 
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were  tested  and  evaluated  for  the  two  types  of  equations  solved  in  this  study,  namely, 
transport  and  Poisson  equations.  For  those  geometries  examined  that  require  an 
outflow  boundary,  a  buffer  domain  technique  was  developed  for  the  vorticity-stream- 
function  formulation  that  allows  disturbances  and  vortex  structures  to  exit  out  of  the 
domain  of  interest  without  reflection. 

The  spectral  method  was  first  tested  using  the  shear-driven  cavity,  both  standard 
and  regularized,  to  validate  the  method,  ascertain  its  accuracy,  and  address  the  issues 
mentioned  above.  Results  for  various  Reynolds  numbers  and  grid  discretization  were 
presented  along  with  favorable  comparisons  between  the  two  formulations  tried  in 
this  study  and  with  other  numerical  results.  The  multidomain  method  was  also 
examined  as  to  its  performance  with  the  shear-driven  cavity  problem.  In  particular, 
the  proper  choice  of  interface  conditions  was  investigated  for  the  two  types  of 
equations  being  used.  Next,  the  simulation  of  the  linear-stability  problem  on  a  flat 
plate  boundary  layer  was  examined  to  test  the  accuracy  of  the  interface  conditions 
with  respect  to  the  passage  of  small  amplitude  waves.  This  problem  provides 
validation  for  the  spectral  multidomain  method  for  use  in  studying  receptivity,  and 
also  introduced  additional  complexities  such  as  grid  stretching,  and  far-fieldfield  and 
outflow  boundary  conditions.  The  spectral  technique  in  generalized  coordinates  is 
given  by  Blodgelt,  K.Ghia  and  Street(1994). 

The  present  method  was  tested  to  insure  that  the  base  flow,  the 
Blasius  boundary-layer  solution,  was  maintained  in  the  absence  of  disturbances.  Two 
cases  of  linear  stability  were  simulated.  The  results  were  compared  with  other  DNS 
results  and  linear  stability  theory.  Grid  studies  were  performed  and  the  results  were 
checked  to  insure  that  they  are  linear.  Two  additional  geometries  were  examined  that 
make  use  of  the  generalized  curvilinear  coordinates  and  for  which  leading-edge 
receptivity  is  investigated:  the  parabola,  and  the  flat  plate  with  an  elliptic  leading 
edge.  Base-flow  results  were  obtained  and  excellent  comparisons  were  observed 
between  the  present  parabola  results  and  those  of  Davis  (1972).  Good  qualitative 
comparisons  were  obtained  between  the  results  from  the  present  method  and  the 
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finite-difference  results  of  Lin  et  al  (1992).  A  receptivity  study  was  attempted  on  the 
parabola  geometry,  but  the  problem,  when  discretized,  proved  to  be  ill-conditioned. 
This  led  to  unacceptable  errors  and  numerical  instabilities.  However,  receptivity  to 
plane  acoustic  waves  was  successfully  simulated  for  the  flat  plate  with  and  elliptic 
leading-edge  geometry.  Favorable  comparisons  with  the  results  of  Lin  et  al.  were 
obtained  for  two  values  of  the  leading-edge  aspect  ratio  and  two  values  of  the  forcing 
frequency.  In  all  of  the  receptivity  results,  the  Tollmien-Schlichting  wave  was  seen  to 
appear  downstream  of  the  ellipse-plate  juncture.  Some  preliminary  results  on 
receptivity  were  discussed  by  Blodgett,  K.Ghia,  U.Ghia  and  Streett(1994). 

The  next  two  studies  by  Duncan  and  Noll  are  currently  still  being 
investigated  and  they  will  now  be  discussed. 


2C  Direct  Numerical  Simulation  of  Transition  and  Turbulence  In 
Complex  Geometry  Internal  Flows  -  Bradley  Duncan,  Ph.D  student 

Purpose 

The  simulation  of  transition  to  turbulence  as  well  as  fully  developed  turbulence  is 
key  to  producing  accurate,  predictive  simulations  of  high  Reynolds-number  flows. 
Correct  simulation  of  turbulent  unsteady  flows  over  aerospace  vehicles  requires  that 
fine-scale  boundary  layer  turbulence  is  correctly  computed,  so  that  the  large-scale 
effects  of  boundary  layer  separation  and  resulting  unsteadiness  will  also  be  accurate. 
In  the  discipline  of  direct  numerical  simulation  (DNS)  of  flows  governed  by  the 
incompressible  Navier-Stokes  equations,  all  of  the  fundamental  scales  of  transitional 
and  turbulent  phenomena  are  computed  in  the  simulation.  The  issues  that  are  then 
faced  are  primarily  computational,  i.e.,  the  computation  is  intensive  and  complex. 
The  addition  of  geometric  complexity  exacerbates  the  computational  difficulties.  In 
order  to  eventually  attain  the  capability  to  perform  simulations  of  turbulent  unsteady 
flow  over  a  three-dimensional  wing,  the  purpose  of  the  current  study  is  to  work 
toward  DNS  methodology  for  complex  geometries. 
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The  present  work  is  focused  on  resolving  computational  and  numerical  analysis 
issues  in  DNS  of  transition  and  turbulence  in  flows  with  some  geometric  complexity, 
but  within  the  domain  of  application  of  traditional  spectral  methods.  Spectral 
methods  are  useful  for  creating  benchmark  results  with  which  other  less  accurate  but 
more  flexible  numerical  approaches  can  be  validated.  Large  eddy  simulation(LES) 
sub-grid-scale  models  are  also  developed  using  spectral  methods,  then  applied  to 
other  approaches.  The  goal  of  the  current  project  is  to  develop  a  general  DNS 
approach  for  spectrally-accurate  simulation  in  complex  geometries,  that  is,  geometries 
with  arbitrary  (but  smooth)  coordinate  transformations. 

The  current  application  of  the  approach  is  to  DNS  of  instability  and  transition  in 
periodic  channel  flows.  The  work  has  focused  on  the  development  of  the 
computational  analysis,  and  has  been  successfully  tested  on  a  curved  channel 
geometry.  The  final  objective  of  the  current  work  is  to  obtain  benchmark  results  of 
transition  to  turbulence  in  the  more  complicated  geometry  of  channel  flow  with 
distributed,  periodic,  surface  roughness.  These  benchmark  results  will  be  of  interest 
to  the  transition  community  as  it  investigates  the  receptivity  of  channel  flow  to  surface 
roughness  and  the  corresponding  route  to  turbulence. 

Progress 

The  study  is  currently  underway.  Significant  progress  had  been  made  in  DNS  of 
instability  and  transition  in  periodic  channel  flows.  The  work  has  focused  on  the 
development  of  the  computational  analysis,  and  has  been  successfully  tested  on  a 
curved  channel  geometry,  a  two-dimensional  channel  with  large-amplitude  waviness, 
and  a  three-dimensional  channel  with  small-amplitude  waviness.  Key  milestones  have 
been  achieved  in  the  present  study  in  the  specific  areas  of  algorithm  development, 
programming  methodology,  and  code  validation.  DNS  results  have  been  obtained  for 
several  cases. 
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Algorithm  development 

The  development  of  spectral  methods  for  complex  geometry  has  been 
challenging,  and  has  mandated  new  strategies  than  have  been  used  in  previous  studies 
with  traditional  spectral  methods.  In  order  to  improve  numerical  stability  and 
robustness,  two  fully-coupled  schemes  were  developed  for  the  study.  Both  schemes 
avoid  the  commonly-used  time-splitting  method  for  decoupling  velocity  and  pressure. 
Instead,  the  velocity  and  pressure  systems  are  solved  within  a  multi-level  iteration 
scheme,  and  the  fully-coupled  linear  algebra  system  is  converged.  The  two  variants 
on  the  fully-coupled  scheme  are  obtained  by  using  either  a  fully-implicit  second-order 
time-discretization  or  a  semi-implicit  second-order  time-discretization  in  which  the 
convective  terms  are  treated  with  an  explicit  Adams-Bashforth  discretization.  The 
fully-implicit  scheme  employs  Newton’s  method  as  the  outer  iteration  level,  and  for 
each  iteration,  a  linearized  system  is  solved.  For  the  semi-implicit  scheme,  iteration  is 
not  required  at  the  outer  level.  For  both  methods,  the  linear  system  is  solved  using 
preconditioned  GMRES  iteration.  Preconditioning  is  performed  using  an 
approximate  lower-upper  (LU)  factorization  of  the  equations  which  decouples  the 
velocity  and  pressure  solution  procedures.  The  velocity  and  pressure  solutions  are 
also  obtained  using  preconditioned  iterative  methods.  The  pressure  system  uses  a  fast 
Poisson  solver  technique  obtained  by  a  transformation  into  wavenumber  space.  The 
DNS  strategy  has  been  significantly  refined  in  the  last  six  months,  and  has  been 
successfully  implemented  into  an  efficient  solution  procedure  for  DNS  of  complex 
geometry  channel  flows. 

Programming  methodology 

In  any  numerical  simulation  code  development  effort,  successful  design  and 
implementation  of  the  details  of  the  programming  approach  are  critical  to  the  success 
of  the  project.  Currently,  research  is  underway  in  the  Computational  Fluid  Dynamics 
Research  Laboratory  to  develop  object-oriented  numerics  (OON)  technology  which 
allows  the  rapid  prototyping  and  development  of  complex  CFD  applications.  The 
fundamental  idea  of  OON  is  the  re-use  of  numerical  algorithms  that  are  universal  in 
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nature  and  are  expressed  in  mathematical  language.  In  the  OON  approach,  these 
algorithms  also  must  be  implemented  in  an  optimized  manner  such  that  computing 
performance  meets  or  exceeds  that  obtained  with  a  conventional  non-object-oriented 
approach.  While  this  effort  is  in  its  early  stages,  the  current  development  of  DNS  of 
channel  flows  is  an  experimental  case  study  for  the  application  of  OON  technology. 
The  code  has  been  developed  in  a  completely  object-oriented  framework,  with 
interchangeable  components  each  representing  a  single  abstraction  of  the  CFD 
problem.  This  allows  much  needed  experimentation  with  different  numerical 
approaches  at  both  the  macro-  and  micro-levels  of  the  code  structure.  Individual 
DNS  “applications”,  or  codes,  can  be  developed  rapidly  from  the  existing  components 
in  order  to  try  a  new  approach.  Within  the  last  six  months,  this  effort  has  been 
successful  in  allowing  testing  of  various  numerical  approaches,  and  in  selection  of  an 
efficient  solution  procedure  for  DNS  of  complex  geometry  channel  flows. 

Code  validation 

During  each  stage  of  development  of  the  current  DNS  code,  detailed 
numerical  validations  have  been  performed.  First,  analytically-prescribed  solutions 
have  been  used,  along  with  a  forcing  term  in  the  equations,  to  test  the  convergence  of 
the  residual  and  computed  solution  with  increasing  spatial  and  temporal  resolution. 
This  approach  will  continue  to  be  used  as  grids  for  new  geometries  are  examined  and 
as  changes  to  the  algorithms  are  introduced,  in  order  to  insure  that  accuracy  has  not 
been  compromised.  Second,  validation  cases  involving  the  evolution  of  linear  waves 
in  plane  and  curved  channel  flows  have  been  used  to  examine  the  accuracy  of  the 
solution.  The  growth  rate  and  frequency  of  oscillation  of  these  waves  compare  well 
with  the  linear  stability  values,  and  allow  detailed  study  of  the  effects  of  spatial  and 
temporal  resolution  and  wave  amplitudes.  Finally,  long-time  computations  of  two- 
dimensional  saturated  states  in  plane  and  curved  channel  flow  have  demonstrated  that 
the  nonlinear  evolution  of  waves  is  computed  correctly.  Saturated  states  have  been 
computed  for  Tollmien-Schlicting  (TS)  waves  and  Dean  vortices  in  channels  with 
different  Reynolds  number  and  curvature.  Most  of  the  validation  stages  were 
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performed  early  in  the  study;  however,  the  validations  are  repeated  as  new  code 
developments  are  applied.  Also,  the  computational  time  required  to  compute  two- 
dimensional  saturated  states  has  recently  been  greatly  reduced  due  to  algorithm 
improvements,  and  results  have  been  obtained,  for  Dean  vortex  flows,  which  compare 
well  to  other  published  results. 

Results 

Though  the  code  is  not  yet  in  “production”  mode,  recently,  three-dimensional 
results  have  been  obtained  which  demonstrate  the  code’s  potential  for  calculation  of 
transition  and  turbulence.  Starting  with  a  two-dimensional  saturated  TS  wave  plus  a 
small  streamwise  vortex  disturbance,  the  evolution  of  a  highly  three-dimensional  flow 
containing  a  single  lambda  vortex  was  obtained  in  a  curved  channel  flow.  Subsequent 
coarse  grid  results  show  the  breakdown  of  large  scale  structure  of  the  flow  into  a 
broadband  small-scale  structure  indicative  of  turbulence.  The  lambda  vortex  structure 
compares  well  with  expected  results,  and  the  subsequent  breakdown  to  turbulence 
demonstrate  correct  behavior  of  the  simulated  results,  though  the  small-scale 
structures  are  not  fully  resolved.  Figure  3  shows  the  lambda-vortex  at  a  time  of  *=75. 
Streamlines  are  shown  lifting  off  from  the  bottom  wall  and  wrapping  around  the 
vortex  (the  head  of  the  vortex  is  on  the  domain  boundary).  On  the  wall  and  the 
inflow  and  side  domain  boundaries,  contours  of  the  magnitude  of  vorticity  are  shown. 
The  wall  vorticity  shows  the  angular  structure  of  the  lambda  vortex,  and  the  spanwise 
non-uniformity  can  be  seen  in  the  inflow  plane.  The  side  domain  boundary 
predominantly  shows  the  characteristic  TS  wave  which  is  undergoing  three- 
dimensional  breakdown. 

Two-dimensional  wavy-channel  results  are  shown  in  Figure  4  as  contours  of 
the  x-component  of  velocity,  the  pressure,  and  the  vorticity  magnitude.  The  results 
are  shown  at  and  early  time  of  *=6.5,  and  were  started  with  an  initial  condition 
obtained  by  projecting  a  fully-saturated  TS  wave  onto  a  divergence-free  solution 
space.  The  large-scale  waviness  adds  a  vortex  shedding  event  to  the  dynamics  of  the 
two-dimensional  instability,  as  can  be  seen  in  the  vorticity  plot. 
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Dean  vortex  solutions  were  obtained  for  two-dimensional  flow  through  the 
channel  cross-section.  At  a  Reynolds  number  of  Re= 250,  this  two-dimensional  flow 
in  unstable  to  streamwise  disturbances.  A  low-amplitude  surface  waviness  was  used 
to  excite  this  instability,  and,  in  Figure  5  results  are  shown  at  t=l,  where  a  streamwise 
variation  in  the  Dean  vortex  rolls  can  be  seen. 

Presentations 

The  code  development  research  as  well  as  some  validation  and  demonstration 
results  have  been  recently  discussed  by  Duncan  and  Ghia(l 997a,  1997b)  at  two 
conferences,  the  First  AFOSR  International  Conference  on  DNS/LES  at  Louisiana 
Tech  University  in  August,  1997,  and  at  the  American  Physical  Society  Division  of 
Fluid  Dynamics  meeting  in  San  Francisco  in  November. 

Proposed  Work 

The  future  goal  of  this  study  are  to  obtain  benchmark  results  for  three- 
dimensional  instability  of  flows  in  curved  channels  with  wall  waviness.  Both  TS-wave 
and  Dean-vortex  simulations  are  currently  being  performed,  and  the  final  results  will 
be  presented  at  conferences,  submitted  for  publication,  and  will  be  the  culmination  of 
Brad  Duncan’s  Ph.D.  research. 


2  D  Development  of  Object  Oriented  Software  For  The  Study  Of 
Compressibility  Effects  on  Dynamic  Stall-  Chris  Noll,  Ph.D  student 


Dynamic  Stall 

One  of  the  goals  for  this  study  was  to  construct  a  simulation  tool  for  studying  the 
compressibility  effects  on  the  dynamic  stall  of  an  airfoil.  Compressibility  effects,  such 
as  shocks,  already  become  prevalent  at  free  stream  Mach  numbers  as  low  as  0.2. 
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Because  of  the  unsteady  flow  induced  by  the  maneuvering  body,  local  Mach  numbers 
can  exceed  unity  in  this  case.  Compressibility  effects  have  the  effect  that  they  lessen 
the  advantages  of  a  maneuver.  As  the  Mach  number  increases,  the  maximum  lift 
coefficient  Ci  decreases,  and,  the  dynamic  stall  occurs  at  lower  angles  of  attack. 
From  the  experimental  studies  of  Chandrasekhara  Carr  and  Wilder(1993  ),  this  was 
attributed  to  two  effects.  At  lower  free  stream  Mach  number  (from  0.2  to  0.4),  a 
shock  forms  above  the  airfoil  surface  triggering  leading-edge  separation,  normally 
associated  with  incompressible  dynamic  stall  at  high  Re.  At  higher  free-stream  Mach 
numbers  (above  0.4),  the  formation  of  numerous  small  shocklets  is  thought  to 
increase  the  local  vorticity  production  hastening  the  dynamic  stall  event. 

Through  initial  simulation  studies,  several  difficulties  were  encountered.  Since 
the  simulation  software  developed  is  time-accurate,  it  was  not  possible  to  impulsively 
start  the  airfoil  from  rest.  When  this  was  done,  a  shock  would  form  enveloping  the 
airfoil  and  then  rapidly  move  out  from  the  airfoil  as  time  progressed.  This  feature 
caused  many  numerical  problems  as  the  code  was  not  yet  designed  to  handle  moving 
shocks.  A  scheme  was  devised  to  accelerate  the  body  from  rest.  This  was  beneficial 
because  the  modifications  to  the  fluid  dynamics  equations  to  allow  for  a  moving  grid 
were  necessary  in  order  to  compute  the  pitching  airfoil  case. 

Another  difficulty  that  was  introduced  as  a  result  was  with  the  traditional 
definitions  of  Mach  number  and  Reynolds  number.  Since  the  body  is  initially  starting 
from  rest  and  remains  moving  at  a  slow  motion  over  the  first  characteristic  time,  the 
Mach  number  and  Reynolds  numbers  of  the  flow  begin  at  zero  and  slowly  move 
toward  their  final  values  progressing  through  the  flow  regimes  in  between. 

Finally,  after  the  airfoil  was  kept  at  a  steady  speed  for  approximately  30 
characteristic  time,  it  was  pitched  up.  Oscillations  developed  in  the  flow  wherever 
flow  features  (such  as  the  wake)  passed  through  regions  where  the  grid  was  clustered. 
This  problem  is  still  unresolved,  but  it  is  believed  that  it  is  due  to  the  absence  of 
artificial  dissipation  in  the  model  equations.  The  model  uses  a  non-traditional 
staggered  continuity  equation.  This  is  being  done  in  order  to  remove  the  extra 
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boundary  condition  on  pressure,  namely  that  the  pressure  gradient  normal  to  the  body 
surface  is  zero,  as  this  condition  does  not  hold  when  the  flow  separates  from  the 
surface. 


Streamlines  and  velocity  vectors,  colored  by  local  Mach  number,  appear  in  the 
figure  below.  This  is  an  initial  result  for  a  coarse-grid  case  and  has  the  oscillation 
problem  discussed  above.  These  initial  results  were  presented  by  Noll  and  K  Ghia 
(1996).  More  recently  a  different  approach  has  been  taken  to  eliminate  the  pressure 
boundary  condition.  The  result  of  this  work  and  a  complete  discussion  of  the 
approach  will  appear  in  Noll  (1998). 

Partial  Differential  Equation  Model 

The  governing  equations  that  are  being  solved  have  been  derived  from  the  first 
principles  of  conservation  of  mass,  momentum  and  energy.  A  time-dependent 
transformation  is  used  to  allow  for  the  motion  of  the  body.  The  Newtonian  form  of 
the  viscous  stress  tensor  is  being  used  along  with  Sutherland’s  law  for  viscosity. 
Stokes’  hypothesis  is  used  for  the  second  coefficient  of  viscosity  and  the  assumption 
of  constant  Prandtl  number  is  used  to  obtain  the  conductivity.  In  addition,  the 
equation  of  state  for  a  perfect  gas  is  used. 

Governing  Equations 


j(/p)t  +  div(pU)  =  0 

j{jpf)x  +  di^pVU  +  {p-  ~KdivV)I  -  \xdefv\  =  0 

7(jpe0)x  +  div[pe0U  +  (p  -  Xdtvf)v  -  \xV  •  dejV  -  kgradT]  =  0 
J=x^-x^ 

where  U  =  V  -W  and  W  is  the  grid  velocity. 
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The  boundary  conditions  that  are  being  used  include  free-stream  flow  far  from 
the  body  and  the  standard  no-slip  and  constant  temperature  condition  on  the  body. 
Much  care  has  been  given,  not  to  introduce  the  commonly  used,  but  inaccurate, 
normal  pressure  gradient  equal  to  zero  condition. 

Numerics 

The  governing  partial  differential  equations  were  discretized  using  central 
differences  for  all  spatial  derivatives  and  a  second-order  trapezoidal  integration 
method  for  the  time  derivative.  This  discretized  system  has  been  solved  using  an  O- 
grid  topology  and  care  has  been  given  to  ensure  that  the  proper  fluid  dynamics 
equations  have  been  solved  on  the  branch  cut. 

The  resulting  non-linear  algebraic  system  is  solved  through  an  approximate 
Newton’s  method.  The  coefficient  matrix  is  loaded  and  frozen  in  time,  so  that  it  can 
be  factored  only  once  for  a  given  number  of  time  steps.  When  convergence  becomes 
difficult,  due  to  the  inaccuracies  thus  introduced  into  the  Jacobian  matrix,  it  is  re¬ 
computed  and  re-factored  at  that  time.  It  was  found  that  the  same  Jacobian  matrix 
could  be  used  over  the  course  of  several  characteristic  time  of  the  computation.  This 
was  true  during  the  constant  angle-of-attack  phase  of  a  computation.  The  Jacobian 
matrix  had  to  be  re-computed  more  often  during  the  pitch-up  phase. 

The  linear  system  resulting  from  the  use  of  Newton’s  method  was  solved  with 
an  approximate  factorization  method  similar  to  Douglas  and  Gunn’s  ADI  method. 
Several  of  these  iterations  were  used  for  each  pseudo-Newton  iteration  until  the 
solution  converged  to  machine-zero  for  each  time  step. 

Object-Oriented  programming 

An  object-oriented  approach  was  taken  to  develop  the  CFD  software.  This 
approach  was  chosen  because  of  the  flexibility  that  it  gives.  This  was  deemed 
necessary  because  many  numerical  methods  were  to  be  tested  and  rapid  development 
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of  the  software  needed  to  solve  the  equations  was  essential  to  allow  for  the  variety  of 
methods  to  be  tested.  The  main  use  for  object-oriented  programming  was  in  the 
solution  of  the  equations.  A  linear  algebra  class  library  was  developed  to  provide 
facilities  for  building  and  solving  multi-dimensional  matrices  of  either  a  dense  or 
structured  form.  The  member  functions  made  extensive  use  of  BLAS  subroutines  to 
get  optimum  performance  for  a  given  hardware  platform.  A  discussion  of  this  work 
was  presented  by  Noll  et  al.  (1996)  and  Brandes-Duncan  et  al.  (1996). 
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Figure  1 


Comparison  of  numerical  results  with  experimental  data  of  Bouard  and 
Coutanceau:  Reo  =  9500,  grid  size  (705  x  91). 

(a)  tD  =  1.0  and  (b)  tD  =  1.5. 
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Wall  Disturbance  Vorticity 
Rat  Plate  With  Elliptic  Leading  Edge,  AR  =  9.0 


Figure  2a  Instantaneous  Disturbance  Vorticity  Along  the  Surface  Re  =  2400 
AR  =  9,  F  =  230x10^. 


Disturbance  Stream  Function 
Rat  Plate  with  Elliptic  Leading  Edge,  AR  =  9.0 


Figure  2b  Instantaneous  Stream  Function  Contour  Plot,  Re  =  2400,  AR  =  9, 
F  =  230  x  10-6 
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Figure  3.  Integrated  streamlines  starting  near  the  lower- wall  of  a  curved  channel,  and  vorticity  magnitude 
contours  on  the  lower  wall,  inflow  plane,  and  side  plane.  Flow  is  from  upper  left  to  lower  right.  The 
parameters  for  the  simulation  are  Re=5000,  radius  of  curvature  (/?c)=100,  a=100,  P=1.25,  grid 
resolution=32x32x64,  Af=0.01,  /=75. 


Figure  4.  Contours  of  y-component  of  velocity  for  a  twisting  Dean  vortex  developing  in  a  curved  channel 
with  small  amplitude  (a=0.01)  surface  waviness.  The  flow  is  from  upper  left  to  lower  right.  The 
parameters  of  the  simulation  are  Re=250,  Rc=19 ,  a=135,  p=2.5,  Af=0.01,  t=  1,  grid  resolution=  16x1 6x32. 


Figure  5.  Contours  of  a)  x-component  of  velocity,  b)  pressure,  and  c)  vorticty  magnitude,  all  for  a  TS  wave 
developing  in  a  wavy  channel.  Flow  is  from  left  to  right.  The  amplitude  of  the  channel  wall  is  0=0.1,  and 
the  parameters  of  the  simulation  are  Re=5000,  Rc=  100,  a=100,  Af=0.01,  t= 6.5,  grid  resolution  =  48x64. 
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Figure  6.  Streamlines  and  velocity  vectors  colored  by  Mach  contour  for  pitching  NACA  0012  at  Mach 
number,  M=0.2,  and  Reynolds  number  1000. 
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Abstract 

A  vorticity-stream  function  circulation 
Too)  formulation  is  explored  as  an  extension  of 

the  earlier  work  of  the  authors  to  focus  upon  the 
accurate  implementation  of  the  far-field  boundary 
condition  for  the  study  of  two-dimensional 
incompressible  unsteady  flows  past  Joukowski  airfoils. 
The  distinction  in  the  implementation  of  the 
formulation  is  achieved  through  the  incorporation  of 
the  viscous  circulation  .  The  resulting  system  of 

governing  equations  is  solved  numerically  using  the 
Block  Gaussian  Elimination  (BGE)  procedure  to  solve 
the  integral  form  of  the  curl  equation,  and  the 
Alternating  Direction  Implicit  (ADI)  technique  of 
Douglas  and  Gunn  to  solve  the  vorticity  transport 
equation,  with  an  iterative  procedure  used  to  solve  the 
coupled  governing  equations  and  boundary  conditions 
as  well  as  ensure  pressure  continuity  along  the  airfoil 
surface.  Due  to  a  lack  of  available  data  pertaining  to 
Joukowski  airfoils,  the  validation  study  is  performed 
using  the  circular  cylinder  geometry,  considering  the 
wake  development  at  early  time  for  two  different  values 
of  the  Reynolds  Re  number.  Additionally,  the  grid 
independence  of  the  solutions  obtained  is  established 
using  the  Grid  Convergence  Index  ( GCJ ).  Additional 
parallel  studies  are  currently  being  conducted  by  the 
authors  at  Computational  fluid  Dynamics  Research 
Laboratory  (CFDRL)  to  examine  the  wall  vorticity 


boundary  condition  and  compare  the  force  coefficient 
data  obtained  with  the  available  experimental  data. 

1.  INTRODUCTION 

The  understanding  of  the  flow  mechanisms 
inherent  in  unsteady  high-ite  flows  past  lifting  airfoils 
lends  valuable  information  towards  the  aerodynamic 
design  of  these  devices.  However,  these  flows  have 
been  found  to  exhibit  a  variety  of  complex  flow 
phenomena.  Much  effort  has  been  expended  by 
researchers  in  attempting  to  investigate  these  flows. 
However,  experiments  involving  high-/te  flows  are,  in 
general,  very  costly  and  are  often  unable  to  provide  the 
detailed  quantitative  information  concerning  the 
temporal  behavior  of  the  flow  field.  Additionally,  even 
with  the  advent  of  large-memory  high-speed 
supercomputers,  the  task  of  numerically  simulating 
high-ite  separated  flows  over  the  complex 
configurations  which  are  of  practical  interest  remains 
formidable,  due  to  the  widely  disparate  length  scales 
present  in  these  flows. 

An  alternative  to  the  prohibitive  costs 
associated  with  high-ite  flow  simulations  lies  in  the 
direct  numerical  simulation  of  low-  and  moderat z-Re 
flows  using  the  Navier-Stokes  equations.  The 
presently-available  supercomputer  resources  have  made 
the  calculation  of  these  flow  solutions  practical,  since 
they  typically  require  considerably  less  memory  and 
CPU  resources  than  their  higher-Re  counterparts. 
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Since  these  lower-ite  systems  often  exhibit  locally  the 
nonlinearities  characteristic  of  high  -Re  separated 
flows,  detailed  examination  of  these  flow  solutions  may 
provide  the  insight  which  is  needed  for  understanding 
of  high-Re  flow  systems.  Towards  this  end,  it  is  the 
purpose  of  the  present  study  to  implement  an 
alternative  far-field  boundary  condition  and  examine 
its  effect  on  the  overall  flow.  To  achieve  this,  a  newly- 
developed  formulation,  presented  in  the  Section  2,  is 
implemented  for  application  to  low-speed  flows. 

Although  the  current  study  focuses  on  the 
development  and  subsequent  validation  of  a  vorticity- 
stream  function-circulation  (#,^,1^)  formulation  for 

application  to  Joukowski  airfoils,  the  ability  of  the 
formulation,  which  uses  generalized  orthogonal 
coordinates,  to  incorporate  an  arbitrary  sequence  of 
conformal  mappings  enables  a  range  of  flow 
configurations  to  be  studied.  The  versatility  of  the 
formulation  makes  it  advantageous  to  consider  a 
validation  study  using  the  circular  cylinder  geometry 
due  to  the  lack  of  available  data  pertaining  to 
Joukowski  airfoils.  This  approach  was  used  previously 
by  Rohling  (1991)  to  validate  the  vorticity-stream 
function  formulation,  which  was  subsequently  applied 
to  conduct  a  detailed  analysis  of  the  flow  past  a 
symmetric  Joukowski  airfoil  at  various  angles  of 
attack;  for  certain  combinations  of  flow  parameters,  the 
flow  exhibited  a  chaotic  behavior. 

The  value  of  the  Reynolds  number.  Re,  has 
been  shown  to  be  of  considerable  importance  in 
determining  the  characteristics  of  the  flow  past  a 
circular  cylinder.  Morkovin  (1964)  and  Marris  (1964) 
have  presented  detailed  reviews  of  the  behavior  of  the 
flow  past  a  circular  cylinder  as  a  function  of  Re.  They 
have  noted  that,  for  subcritical  values  of  Re,  i.e., 
Re  <  100,000  -  130,000,  the  flow  is  characterized 
by  a  laminar  near  wake  consisting  of  coherent  vortex 
structures  which  remain  symmetric  for  a  period  of  time 
and  subsequently  shed  alternately  from  the  rear  of  the 
cylinder.  This  periodic  vortex  shedding  phenomenon 
has  also  been  noted  in  the  computational  investigations 
of  Liu  (1987),  K.  Ghia  et.  al.  (1987)  and  Blodgett 
(1990),  but  at  very  low  Re  values.  Additionally  Braza, 
Chassaing  and  Ha  Minh  (1986)  have  noted  the 
presence  of  periodic  vortex  shedding  for  relatively  low 
values  of  Re  and  believe  this  to  be  an  intrinsic  flow 
phenomenon  which  is  contained  in  the  Navier-Stokes 
equations. 

Of  particular  relevance  to  the  current  study  is 
the  experimental  investigation  of  Bouard  and 


Coutanceau  (1980),  which  considered  the  early  wake 
development  of  the  flow  past  an  impulsively-started 
circular  cylinder  at  low-,  moderate-  and  high-Reynolds 
numbers.  Their  study  utilized  a  flow  visualization 
technique,  which  permits  detailed  observation  of  the 
temporal  behavior  of  the  flow  structures,  to 
successfully  examine  both  quantitatively  and 
qualitatively,  the  evolving  flow  structures  and  vortex 
interactions.  Their  investigations  of  the  flow  at 
moderate  Re  demonstrated  a  strong  qualitative 
agreement  with  the  results  of  Honji  and  Taneda  (1969), 
Taneda  (1972)  and  Thoman  and  Szewczyk  (1969). 
Quantitatively,  their  consideration  of  separation  bubble 
lengths  at  early  times  were  in  close  agreement  with  the 
results  of  Collins  and  Dennis  (1973)  and  Panikker  and 
Lavan  (1975),  among  others.  It  is  therefore  chosen  to 
utilize  the  experimental  results  of  Bouard  and 
Coutanceau  (1980)  to  validate  the  results  of  the  present 
study. 

2.  Mathematical  Formulation 


2.1  Governing  Equations 


The  earlier  analysis  of  Osswald,  K.  Ghia  and 
U.  Ghia  (1985),  which  considered  the  numerical 
solution  of  the  Navier-Stokes  (NS)  equations  in  terms 
of  the  derived  variables  of  vorticity  co  and  stream 
function  y/ ,  has  been  modified  to  incorporate  an 
alternative  far-field  boundary  condition,  which  results 
in  the  coupling  of  the  viscous  circulation  with  the 

disturbance  stream  function  y/D .  Additionally,  in  an 

on-going  study  at  CFDRL,  Osswald  has  investigated 
the  advantages  of  obtaining  the  stream  function  ys 

from  an  integral,  rather  than  differential,  relation, 
which  facilitates  proper  treatment  of  the  weak 
singularity  existing  at  the  body  surface.  Then,  the 
governing  equations  take  the  form  of  the  temporally- 
parabolic  vorticity  transport  equation  and  the  spatially- 
elliptic  integral  form  of  the  curl  equation  for  velocity, 
also  known  as  Stokes  Theorem.  These  are  given  in 
terms  of  generalized  orthogonal  coordinates  as: 

Vorticity  Transport  Equation: 
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Integral  form  of  Curl  Equation: 

ii^a,dM2=i(-j=-%n-^^re2\clI  (2) 
dA  e\'i8  d?  Jga*1  ) 

with  V  •  V  =  0  being  automatically  satisfied  by  the 
stream  function  formulation. 

Since  the  stream  function  does  not  appear 
explicitly  in  Eq.  (1),  it  is  necessary  to  incorporate  the 
relation  between  stream  function  and  velocity  into  the 
formulation,  given  as 

V=Vyxk  ,  (3) 

where  k  denotes  the  unit  vector  normal  to  the  two- 
dimensional  flow  domain.  Osswald  has  show  that 
when  Eqs.  (2-3)  are  discretized  using  a  conservative 
central-difference  scheme,  the  resulting  set  of  algebraic 
equations  is,  for  points  which  lie  in  the  interior  of  the 
flow  domain,  exactly  the  same  as  that  obtained 
previously  by  Osswald  et  al.  (1985),  while  having 
considered  the  discretized  form  of  the  spatially-elliptic 
stream  function  equation,  given  as 

vV  =  -0)  .  (4) 

2.2  Boundary  and  Initial  Conditions 

The  numerical  solution  of  the  temporally- 
parabolic  vorticity  transport  equation  and  the  spatially- 
elliptic  integral  form  of  the  curl  equation  necessitates 
the  specification  of  appropriate  boundary  conditions, 
both  at  the  far-field  boundaries  of  the  flow  domain  and 
at  the  airfoil  surface,  as  well  as  an  appropriate  initial 
condition.  Previous  experience  with  maneuvering 
simulations  has  revealed  the  critical  role  of  the 
boundary  conditions  in  determining  the  flow  evolution. 
The  present  formulation  represents  an  exploratory 
study  where  boundary  conditions  are  re-examined  to 
achieve  continuity  in  the  surface  pressure  distribution. 

Previously,  Rohling  (1991),  in  his  analysis  of 
the  unsteady  incompressible  flow  past  a  12%  thick 
symmetric  Joukowski  airfoil  using  the  {(o,y/) 

formulation,  implemented  boundary  conditions  which 
consisted  of  the  no-slip  condition  along  the  body 
surface  and  free-stream  conditions  at  the  far-field 
boundaries.  However,  recently,  Osswald,  while 
developing  an  analysis  to  study  the  three-dimensional 
flow  past  a  finite  wing  using  the  incompressible 
unsteady  NS  equations,  found  the  Cauchy-Riemann 


operator  to  be  singular  for  multiply  connected  domains, 
regardless  of  whether  a  C-grid  or  an  O-grid  is 
employed.  This  results  in  a  single  degree  of 
degeneracy  which  needs  to  be  accounted  for.  Since  the 
no-slip  condition,  applied  at  the  airfoil  surface,  is  well- 
behaved,  and,  since  the  vorticity  at  infinity  is 
identically  zero,  and  hence,  clearly  well-behaved,  the 
asymptotic  behavior  of  the  far-field  stream  function 
was  carefully  examined  in  hopes  of  resolving  the 
singularity  in  the  Cauchy-Riemann  operator.  The 
study  of  Yang  (1992),  which  included  an  asymptotic 
analysis  of  the  far-field  behavior  of  the  stream  function 
obtained  using  a  Laurent  series  expansion,  revealed 
that  the  viscous  circulation  exhibits  a  similar  behavior 
at  infinity.  In  light  of  this  behavior,  the  original 
(<o,^)  formulation  has  been  modified  to  include  an 

alternative  far-field  boundary  condition  for  the  stream 
function,  which  results  in  the  coupling  of  the  viscous 
circulation  at  infinity  with  the  vorticity  transport  and 
integral  form  of  the  curl  equations.  The  resulting 
(^.Too)  formulation  necessitates  that  additional 

conditions  be  used  to  (/)  provide  the  value  of  roo(f)  at 

each  time  step,  and  (if)  uniquely  specify  its  relationship 
to  the  disturbance  stream  function  at  infinity. 

To  remove  the  numerical  difficulties 
associated  with  the  unbounded  behavior  of  the  stream 
function  at  infinity,  the  stream  function  is  decomposed 
as: 

^  =  ,  (5) 

where  denotes  the  stream  function  corresponding 
to  the  zero-lift,  i.e.,  zero-circulation  inviscid  flow 
solution,  y/G  represents  the  stream  function 
corresponding  to  a  unit  circulation  potential  slip 
solution,  with  the  viscous  circulation  Too  which 

remains  nonzero  at  infinity,  and  y/D  represents  the 

disturbance  stream  function,  which  is  chosen  as  the 
dependent  variable  in  the  numerical  solution  of  the 
integral  form  of  the  curl  equation.  The  use  of  the  zero- 
lift  inviscid  stream  function  y/° ,  rather  than  the 

inviscid  stream  function  y/1 ,  which,  in  general, 

possesses  a  non-zero  circulation,  serves  to  account  for 
the  inviscid  behavior  except  for  that  due  to  the 
circulation  in  the  decomposition  of  the  stream  function. 
It  should  be  noted  that,  since  the  unbounded  behavior 
at  infinity  associated  with  the  zero-lift  inviscid  stream 
function  and  the  viscous  circulation  have  been  removed 
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from  yD ,  the  numerical  solution  of  the  disturbance 
stream  function  constitutes  a  well-posed  problem,  since 
y/D  and  its  boundary  conditions  are  nonsingular  over 
the  entire  flow  domain. 

Therefore,  boundary  conditions  are  needed  on 
y/D .  The  examination  of  the  far-field  behavior  of  the 
stream  function  by  Osswald  revealed  that  the 
disturbance  stream  function,  <yD ,  does  not,  as  was 
previously  believed,  approach  zero  at  infinity,  but 
rather,  it  approaches  an  unknown  constant  Therefore, 
the  Dirichlet  boundary  condition  on  the  disturbance 
stream  function  used  in  the  study  of  Yang  (1992), 
namely,  the  enforcement  of  y/D  =  0  at  infinity,  was 
modified  in  the  current  study  to  reflect  this.  Then,  the 
boundary  conditions  employed  in  the  present  study  are 
given  as: 

At  the  airfoil  surface: 

V  =  0  (No-slip  condition) ,  (6) 

At  infinity: 

\ftD  =  y/^  (Dirichlet  Boundary  Condition) ,  (7) 

where  ^oo  denotes  the  constant  value  of  the  stream 

function  at  infinity.  These  boundary  conditions  are 
shown  graphically  in  Fig.  1.  The  no-slip  condition  is 
implemented  via  the  integral  form  of  the  curl  equation 
applied  over  the  domain  consisting  of  the  half  cell 
surrounding  the  body. 

Previously,  in  the  study  of  Rohling  (1991),  the 
inviscid  potential  flow  solution  at  the  flow  angle  of 
attack  a  f  was  considered  as  the  initial  condition. 

However,  it  was  later  discovered  by  Osswald  that,  for 
the  formulation  considered  in  the  present 

study,  the  use  of  an  inviscid  solution  as  the  initial 
condition  fails  to  provide  the  correct  initial  value  of  the 
viscous  circulation,  and  hence  cannot  be  applied. 
Therefore,  the  initial  condition  considered  in  the 
present  study  is  the  viscous  flow  solution  at  time 
/  =  0+ ,  corresponding  to  an  impulsive  start  from  rest 
of  the  airfoil. 

The  system  of  equations  consisting  of  Eqs.  (1) 
and  (2),  combined  with  their  respective  boundary  and 
initial  conditions,  cannot  be  solved  uniquely  unless  the 


value  of  is  given  at  each  time  step.  As  mentioned 
previously,  specification  of  requires  knowledge  of 

the  viscous  circulation  remaining  at  infinity. 
Therefore,  an  additional  condition  is  needed  which 
gives  the  appropriate  value  of  r<jo(/)  at  each  time  step. 
Mathematically,  this  condition  may  be  written  as 

roo  =  rB-JJ<a  dA  ,  (8) 

D 

where  rB  denotes  the  circulation  at  the  body  surface 
and  D  denotes  the  entire  flow  domain.  Note  that  the 
boundary  condition  <a  =  0  at  infinity  makes  it  possible 
to  impose  the  far-field  limit  of  integration  on  the 
\\<odA  term  at  the  half-cell  location,  rather  than  at  the 

D 

actual  boundary  of  D,  which  lies  at  true  infinity.  The 
difficulties  associated  with  numerical  integration  over 
the  semi-infinite  cells  at  the  far-field  boundary  are  thus 
avoided.  Also,  with  regards  to  the  present  study,  the 
relation 

rB(f)  =  0  ,  (9) 

holds  for  all  time,  since  (/)  the  angular  acceleration  of 
the  airfoil  is  identically  zero,  and  (//)  no  suction  or 
injection  is  applied.  While  Eq.  (8)  provides  the  value 
of  roo(r)  at  each  time  step,  its  relationship  to  the  far- 

field  stream  function  must  be  specified  in  order  to 
ensure  a  well-posed  problem.  This  is  achieved  by 


Figure  1  Boundary  Conditions  for  the  Joukowski 
Airfoil 
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invoking  the  integral  form  of  the  curl  equation  along 
the  body  surface  as 

iV-di^-Tn-WcodA  ,  (10) 

S  B 

where  the  subscript  B  denotes  the  path  one  half  cell 
off  of  the  body  surface  and  B  denotes  the  domain 
between  the  body  and  the  path  denoted  by  B .  The 
above  Eq.  (10),  when  discretized  is  a  linear 
combination  of  the  discretized  form  of  the  no-slip 
conditions  and  thus  there  is  still  a  degree  of  freedom  to 
be  accounted  for.  Following  the  method  of  Yang 
(1992),  this  difficulty  is  resolved  by  eliminating  one  of 
the  discrete  no-slip  conditions  and  imposing  the 
condition  that  the  pressure  be  single-valued  and 
continuous  along  the  airfoil  surface.  Mathematically, 
this  condition  is  given  as 

\dp  =  0,  (11) 

B 

where  the  subscript  B  denotes  the  body  surface,  and 
represents  the  Kutta  condition  for  unsteady  viscous 
flow,  i.e.,  Eq.  (11)  necessitates  that  the  pressure  at  the 
infinitely-sharp  trailing-edge  of  the  airfoil  must  be 
single-valued  and  continuous. 

Equation  (11)  may  be  evaluated 

mathematically  by  considering  the  linear  momentum 
equation  as 

V?  =  Kx5-±(VxS)-^-v(^)  (12) 

and  recalling  that,  at  the  body  surface,  the  boundary 
condition  V  =  0  must  be  satisfied.  Then,  at  the  airfoil 
surface,  Eq.  (12)  is  reduced  to  a  statement  that  the 
variation  of  pressure  along  the  body  surface  depends 
only  upon  the  normal  variation  of  the  vorticity. 
Mathematically,  this  statement  is  given  as 

ldp  =  -l~dl  =  0  ,  (13) 

B  B 

where  n  denotes  the  direction  normal  to  the  body 
surface,  and  dl  is  taken  along  the  body  surface.  In  the 
solution  procedure,  Eq.  (13)  is  combined  with  the  no¬ 
slip  condition  to  provide  a  unique  distribution  of  the 
body  surface  vorticity,  with  pressure  closure  enforced 
at  the  leading  stagnation  point  (LSP).  An  advantage  of 
enforcing  pressure  closure  at  the  LSP  is  that  this  is  a 


point  of  symmetry  for  the  cases  of  (i)  flow  past  a 
symmetric  airfoil  at  zero  angle  of  attack,  and  (ii)  flow 
over  a  non-rotating  circular  cylinder.  Therefore,  for 
these  symmetric  cases,  the  enforcement  of  pressure 
closure  at  the  LSP  results  in  the  inherent  flow 
symmetry  being  retained  throughout  the  calculation. 

3.  Grid  Generation  Procedure 

The  present  study  utilizes  the  grid  generation 
procedure  of  Osswald  et  al.  (1985),  which  uses  a  series 
of  analytical  transformations  to  map  a  unit  square 
computational  domain  into  a  C-grid  topology  in  the 
physical  plane,  with  cubic  spline  transformations  [see 
Rohling  (1991)]  to  provide  precise  grid  point  control. 
This  approach  has  been  shown  by  Rohling  (1991)  to 
result  in  a  robust  and  efficient  grid  generation 
procedure  which  is  capable  of  yielding  nearly  optimal 
grids  for  Joukowski  airfoils  at  various  angles  of 
incidence;  details  of  this  grid  generation  procedure 
may  be  found  in  Rohling  (199 1). 

The  grid  generation  procedure  uses  a  total  of 
four  separate  transformations,  as  shown  in  Fig.  2,  to 
achieve  the  desired  mapping  between  the  physical 
plane  and  the  computational  plane.  The  first  two 
transformations,  which  constitute  a  mapping  of  the 
physical  plane  (Fig.  2  (a))  to  the  complex-potential 
plane  (Fig.  2  (c)),  are  conformal  and  are  representative 
of  the  inviscid  flow  solution,  which,  for  Joukowski 
airfoils,  is  known  in  analytical  form.  Subsequently,  a 
parabolic  conformal  mapping  is  used  to  HunfoldH  the 
complex-potential  plane  about  the  leading-edge 
stagnation  point,  thus  enabling  a  C-grid  topology  to  be 
generated  in  the  physical  plane.  The  final 
transformation  is  merely  orthogonal,  i.e., 
nonconformal,  and  uses  a  contraction  mapping  in  the 
far  field  to  transform  the  doubly-infinite  77-plane  (Fig. 
2  (d))  into  a  unit  square  computational  domain  (Fig.  2 
(e)),  which  makes  it  possible  to  impose  the  far-field 
boundary  conditions  at  true  infinity,  rather  than  at 
some  finite  distance  away  from  the  body.  This 
transformation  also  uses  one-dimensional  cubic  spline 
functions  of  Rohling  (1991)  in  the  interior  region  to 
provide  the  necessary  control  over  the  grid  point 
distribution  by  allowing  the  number  of  cells  in 
specified  stream  wise  and  normal  zones  as  well  as  the 
positions  of  the  zonal  boundaries  to  be  completely 
arbitrary  and  specified  by  the  analyst. 
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(a)  Physical  Plane 


4.  Numerical  Procedure 


(c)  Complex  Potential  Plane 


pL(D  j^1 


uoo 

“LSP1 


LE 


S 


(d)  T[  -Plane 
Uool 


TEU  Branch 


Cut 


n2‘ 

i  =  U  +  i*1 

Lower  TEL 

LSP  TE  Upper 

Branch 

(e)  Computational  Plane 


(0,1) 


(0,0) 


.  U^j 

(1,1) 

II 

UJ' 

TEL  LSP  TEU 

(i.Q) 

The  numerical  solution  procedure  utilized  in 
the  present  study  follow's  the  original  folly-implicit 
time-marching  solution  method  for  unsteady  flows 
developed  by  Osswald,  K.  Ghia  and  U.  Ghia  (1985). 
Modifications  to  this  method  have  been  made  in  order 
to  incorporate  the  recent  advances  of  Yang  (1992)  and 
unpublished  work  of  Osswald  towards  the  correct 
implementation  of  the  boundary  and  initial  conditions. 
A  detailed  description  of  the  solution  procedure  is 
contained  in  Beltz  (1994). 

Although  the  set  of  governing  equations 
presented  in  Section  2  clearly  demonstrates  the 
coupling  of  the  vorticity  and  stream  function  fields,  the 
present  (a^Too)  formulation  contains  a  natural 

separation  of  the  spin  dynamics  of  a  fluid  particle, 
governed  by  the  vorticity  transport  equation,  and  its 
translational  kinematics,  governed  by  the  integral  form 
of  the  curl  equation.  This  permits  the  vorticity  and 
stream  function  fields  to  be  solved  separately. 
Additionally,  the  disturbance  stream  function,  viscous 
circulation  and  surface  vorticity  are  predicted  and 
corrected  via  an  iteration  procedure,  which  involves 
updating  the  surface  vorticity  via  the  no-slip  condition 
and  the  pressure  closure  constraint  This  procedure  is 
repeated  until  convergence  of  the  surface  vorticity  is 
obtained  to  within  a  specified  tolerance.  The  use  of  the 
iterative  procedure  in  obtaining  the  surface  vorticity  at 
the  n  + 1  time  level  results  in  a  formal  accuracy  of  the 
overall  numerical  solution  procedure  of 


In  providing  the  discretization  of  the  vorticity 
transport  equation,  the  original  procedure  of  Osswald, 
K.  Ghia  and  U.  Ghia  (1985)  considered  the  use  of  the 
second-order  accurate  central  difference  scheme  (CDS) 
to  approximate  all  spatial  derivatives.  However,  it  was 
later  shown  by  Yang  (1992)  that  the  flow  solution 
obtained  using  this  differencing  technique  showed 
spurious  oscillations  in  the  vorticity  field  for  the 
Re  =  10,000  flow  past  a  NACA  0015  airfoil 
undergoing  a  constant  pitch-up  motion.  This  was 
attributed  by  Yang  (1992)  to  a  shortcoming  of  the 
central  differencing  technique,  which  retains  second- 
order  accuracy'  but  is  limited  by  the  cell  Reynolds 
number  criterion. 

The  effectiveness  of  artificial  dissipation  in 
eliminating  the  numerically-induced  oscillations 
associated  with  the  CDS  has  been  well-established  by 


Figure  2  Transformations  employed  in  the  grid 
generation  procedure. 
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researchers.  However,  the  amount  of  artificial 
dissipation  required  differs  for  each  specific  case, 
making  its  application  difficult  An  alternative  to  this 
approach  is  the  use  of  a  third-order  accurate  biased 
upwind  difference  scheme  (UDS),  which  has  a  natural 
tendency  to  provide  dissipation,  to  approximate  the 
convective  terms  in  the  vorticity  transport  equation. 
This  approach  was  demonstrated  by  Yang  (1992)  to 
successfully  eliminate  the  spurious  oscillations 
previously  observed,  resulting  in  a  more  accurate 
physical  solution.  Additionally,  Blodgett,  K.  Ghia, 
Osswald  and  U.  Ghia  (1990)  successfully  obtained 
results  for  the  high-/te  flow  past  an  elliptic  cylinder 
using  the  UDS.  Therefore,  in  the  present  study,  the 
diffusion  terms  appearing  in  the  vorticity  transport 
equation  are  discretized  using  a  second-order  accurate 
central  difference  scheme  (CDS),  while,  for  the 
convective  terms,  this  CDS  is  applied  only  in  the  low- 
speed  regions  of  the  flow.  When  the  apparent  flow 
velocity  is  relatively  large,  a  third-order  biased  UDS  is 
applied. 

The  set  of  algebraic  equations  resulting  from 
the  discretization  of  the  vorticity  transport  equation  is 
solved  using  the  alternating-direction  implicit  (ADI) 
technique  of  Douglas  and  Gunn  (1964).  This  method 
possesses  the  advantage  of  computational  efficiency 
and  facilitates  easy  extension  to  the  three-dimensional 
case  of  a  finite  Joukowski  wing.  As  noted  in  Douglas 
and  Gunn  (1964),  application  of  the  ADI  technique 
requires  the  solution  of  two  implicit  pentadiagonal 
matrix  equations,  one  corresponding  to  each 
directional  sweep,  which  are  solved  sequentially  over 
the  entire  flow  domain.  Following  the  method  of  Yang 
(1992),  these  matrix  equations  are  solved  using  the 
generalized  Thomas  algorithm.  The  stability  of  the 
ADI  method  is  assured  provided  that  the 
aforementioned  matrices  are  diagonally  dominant. 
The  use  of  directional  switching  in  the  UDS  serves  to 
consistently  augment  the  diagonal  terms  of  the 
matrices,  thus  increasing  the  stability  envelope  of  the 
ADI  procedure. 

The  linearity  and  constant  coefficients  of  the 
integral  form  of  the  curl  equation  enable  it  to  be  solved 
directly  via  the  Block  Gaussian  Elimination  (BGE) 
technique,  with  the  value  of  the  viscous  circulation 
Too  used  to  determine  the  far-field  stream  function 

.  The  matrix  vector  equation  thus  obtained 
contains  the  block  tridiagonal  matrix  previously 
obtained  by  Osswald  et  al.  (1985),  with  an  extra  row 
and  column  added  to  this  matrix  to  accommodate  the 
coupling  of  the  integral  form  of  the  curl  equation  with 


the  viscous  circulation.  As  shown  by  Beltz  (1994),  the 
presence  of  this  extra  row  and  column  in  the  matrix 
does  not  significantly  alter  the  BGE  solution  procedure 
or  its  computational  efficiency.  The  interior  of  the 
matrix  is  reduced  to  upper-diagonal  form  using  the 
Block  Gaussian  Elimination  (BGE)  procedure  of 
Osswald  et  al.  (1985);  subsequently,  the  techniques  of 
this  method  are  used  to  perform  the  necessary 
elimination  of  the  last  row  and  column.  Since  the 
matrix  consists  entirely  of  constant  metric  coefficients, 
it  is  necessary  to  perform  the  BGE  procedure  only 
once.  A  computationally-efficient  back-substitution 
procedure  is  used  to  compute  y/D  over  the  entire  flow 
domain  at  each  level  of  iteration. 

5.  RESULTS  AND  DISCUSSION 

The  mathematical  formulation  and  numerical 
solution  procedure  presented  here  is  validated  using 
flows  past  the  circular  cylinder  geometry.  This  is 
easily  accomplished  by  eliminating  the  final  Joukowski 
transformation  from  the  sequence  of  analytical 
mappings  used  in  the  grid  generation  procedure,  which 
results  in  the  transformation  of  the  unit  square 
computational  domain  (Fig.  2(e))  to  a  C-grid  topology 
in  the  circular  cylinder  plane  (Fig.  2(b)).  This  minor 
modification  affects  only  the  analytically-calculated 
metric  terms,  and  hence  does  not  impact  the  numerical 
solution  procedure,  thus  allowing  easy  extension  of  the 
current  formulation  to  study  the  unsteady  viscous  flow 
over  a  circular  cylinder. 

As  mentioned  earlier,  the  numerical 
simulation  results  are  validated  using  the  experimental 
study  of  Bouard  and  Coutanceau  (1980),  which 
considered  the  early  stages  of  the  wake  development 
behind  the  cylinder  using  a  flow  visualization 
technique  fine  enough  to  allow  detailed  observation  of 
the  flow  structures.  A  subsequent  analysis  of  the 
photographs  obtained  was  used  to  provide  quantitative 
measurements  of  the  wake  properties,  as  well  as 
qualitative  characteristics  regarding  the  time 
development  of  the  recirculation  region  which  forms 
behind  the  cylinder.  Although  the  geometry  is  simple, 
complex  phenomena  have  been  found  to  occur  at 
higher  Reynolds  numbers. 

The  present  study  considers  flow  past  the 
cylinder  at  two  different  Reynolds  number,  namely,  the 
simple  test  case  of  ReD  -  200  and  also  the  more 
complex  flow  case  of  ReD  =  9500,  where  ReD  denotes 
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the  Reynolds  number  based  on  the  cylinder  diameter, 
defined  as 


2a^oo 

v 


(14) 


where  U ^  is  the  free-stream  velocity,  a  is  the  radius  of 
the  cylinder  and  v  is  the  kinematic  free-stream 
viscosity.  Also,  the  nondimensionalized  time  based  on 
cylinder  diameter  tD  is  taken  as 


(15) 


5,1  Results  of  £^=200  Case 

The  results  of  the  simple  test  case  of 
ReD  =  200  are  presented  in  detail  in  Beltz  (1994),  and 
hence,  only  a  synopsis  of  these  results  is  contained 
here.  The  grid  size  considered  for  the  ReD  =  200  case 
is  (97  x  30) ,  corresponding  to  97  points  in  the 
streamwise  direction  and  30  points  in  the  normal 
direction,  and  a  time  step  of  AtD  =  0.001  is  employed. 
The  use  of  this  relatively  coarse  grid  is  consistent  with 
the  low  Reynolds  number  of  the  flow,  and  the  results 
for  this  case  have  shown  its  ability  to  yield  the  desired 
degree  of  flow  resolution. 


of  the  maximum  wake  width  xw  ID  was  made  and 

compared  with  the  experimental  results.  These 
quantities  are  shown  schematically  in  Fig.  3.  A 
graphical  comparison  of  the  computed  and 
experimental  measurements  of  lw! Dy  w^/D,  and 

xWmm  ID  for  this  case  is  given  in  Fig.  4,  and  these 

results  demonstrate  a  close  agreement  with  the 
experimental  results  of  Bouard  and  Coutanceau  (1980). 
Discrepancies  between  the  numerical  and  experimental 
results  may  be  attributed  to  differences  in  the 
numerical  and  experimental  techniques  used  to 
pinpoint  the  location  of  the  stagnation  streamline. 
Specifically,  experimental  investigations  allow  this 
contour  value  to  be  determined  from  the  flow 
visualization  and  velocity  measurements,  whereas 
numerical  studies  rely  upon  an  interpolation  of  the 
stream  function  values  at  neighboring  grid  points. 
Additionally,  Bouard  and  Coutanceau  (1980) 
acknowledge  that,  despite  the  high  degree  of  flow 
resolution  attained  in  their  study,  a  margin  of  error  as 
great  as  10%  may  be  present  in  their  results  for  the 
measured  wake  properties.  Therefore,  the 
disagreements  observed  between  the  experimental  and 
numerical  measurements  is  considered  to  lie  well 
within  the  realm  of  the  uncertainty  associated  with 
these  data  sets. 

5,2  Results  of  lteI)=9500  Case 


Examination  of  the  vorticity  and  stream 
function  fields  for  this  case  revealed  the  presence  of  a 
distinct  separation  bubble  behind  the  cylinder  at 
tD  =  05  9  indicating  that  flow  separation  has  occurred 

prior  to  this  time.  As  time  progresses,  the  vorticity 
generated  at  the  cylinder  surface  convects  downstream 
into  the  wake,  inducing  a  recirculating  flow  region 
behind  the  cylinder  which  grows  steadily  with  time. 
This  recirculation  region  consists  of  a  single  pair  of 
counter-rotating  vortices  which  remain  symmetric  and 
attached  to  the  cylinder  through  tD  =  3.0 .  At  this  low 
Reynolds  number,  the  wake  development  has  been 
shown  to  occur  without  any  visible  distortion  or 
instability. 

The  above  description  of  the  flow  evolution 
closely  matches  the  qualitative  description  of  the  wake 
development  given  by  Bouard  and  Coutanceau  (1980). 
A  quantitative  comparison  of  (/)  the 
nondimensionalized  wake  length  lw!D,  (//)  the 
nondimensionalized  maximum  wake  width  /  D , 
and  (///')  the  nondimensionalized  streamwise  location 


The  successful  validation  of  the  numerical 
simulation  results  for  the  simple  test  case  of 
ReD  =  200  led  to  the  consideration  of  the  more 

complex  case  of  ReD  =  9500.  A  grid  size  of 
(705x91),  corresponding  to  the  higher  Reynolds 
number  of  the  flow,  was  utilized  for  this  case,  and  the 
time  step  considered  was  tD  =  0.0001 .  The 
corresponding  grid  distribution  is  shown  in  Fig.  5,  and 
the  time  evolution  of  the  vorticity  and  stream  function 
fields  is  shown  in  Fig.  6.  While  the  wake  development 
is  similar  to  the  ReD  -  200  case  in  that  the  flow  is 
mainly  characterized  by  a  symmetric  pair  of  counter¬ 
rotating  vortices  which  form  behind  the  cylinder,  the 
flow  at  ReD  =  9500  is  considerably  more  complex.  As 
shown  in  Fig.  6,  a  secondary  phenomenon  has  clearly 
emerged  prior  to  tD  =  1.0  in  the  stream  function  field 
near  the  wall  at  about  three-quarters  of  the  distance 
between  the  rear  stagnation  point  and  the  separation 
point.  Initially,  this  secondary  phenomenon  consists  of 
an  isolated  secondary  eddy,  which  has  a  rotation 
opposite  to  that  of  the  main  eddy.  In  the  study  of 
Bouard  and  Coutanceau  (1980),  the  formation  of  the 
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secondaiy  eddy  is  attributed  to  the  separation  of  the 
back  flow  from  the  cylinder  wall,  which  is  a  direct 
result  of  the  high  flow  velocity.  As  the  secondary  eddy 
develops  and  grows  in  size,  it  touches  the  boundary  of 
the  main  recirculating  zone  and  splits  into  two  parts, 
forming  another  secondary  eddy  which  is  equivalent  in 
size  and  strength  to  the  first.  This  description  of  the 
flow  evolution  closely  matches  the  behavior  described 
in  the  study  of  Bouard  and  Coutanceau  (1980).  A 
qualitative  comparison  of  the  numerical  and 
experimental  visualization  results,  which  demonstrates 
the  secondary  phenomenon  at  tD  =  1.0  and  tD  =  L5,  is 
given  in  Fig.  7.  This  comparison  clearly  demonstrates 
the  ability  of  the  numerical  simulation  to  accurately 
capture  the  secondaiy  flow  structures  observed  in  the 
experimental  study. 

A  quantitative  comparison  of  the  flow  at 
ReD  =  9500  is  given  in  terms  of  (/)  the 
nondimensionalized  wake  length  lw  l  D  and  (//)  the 
nondimensionalized  circle  plane  coordinates  of  the 
main  eddy  core  / D.b^  /  2D) .  A  comparison  of 

the  numerical  and  experimental  values  of  these 
quantities  given  in  Tables  1  through  3;  the 
corresponding  graphical  comparison  is  given  in  Fig.  8. 
As  for  the  ReD  =  200  case,  discrepancies  between  the 
numerical  and  experimental  results  may  be  attributed 
to  differences  in  the  techniques  used  to  pinpoint  the 
location  of  the  stagnation  streamline.  Note  that  the 
discrepancy  between  the  numerical  and  experimental 
measurements  of  lw!D  increases  with  time,  which 

may  be  expected  due  to  the  increasing  complexity  of 
the  flow.  The  numerical  measurements  of  the 
nondimensionalized  r1 -coordinate  of  the  main  eddy 
core  /  D  are,  in  general,  in  very  close  agreement 
with  the  experimental  results.  Although  the  agreement 
is  not  as  close  for  the  nondimensionalized  r 2  - 
coordinate  /  2D ,  it  should  be  noted  that  the  trend 
is  somewhat  duplicated,  as  shown  in  Fig.  8(c).  In 
general,  the  quantitative  agreement  between  the 
numerical  and  experimental  measurements  of  the  wake 
properties  is  satisfactory;  this  agreement  is  largely  the 
result  of  the  accurate  resolution  of  the  secondary  eddy 
phenomenon  realized  in  the  current  study. 

5.3  Grid  Independence  Study 

The  validation  study  described  thus  far  is  not 
complete  without  an  investigation  of  the  grid 
independence  of  the  solutions  obtained.  This  section 
serves  to  resolve  this  issue  by  analyzing  the  numerical 


simulation  results  obtained  at  ReD  =  9500  using  three 
different  grids.  In  this  analysis,  the  (705x91)  grid 
considered  in  the  previous  section  is  taken  to  be  the 
"medium",  or  "ba seH,  grid  A  coarse  (353  x  46)  grid  is 
obtained  via  a  50%  reduction  of  the  number  of  grid 
points  in  both  the  stream  wise  and  normal  directions 

Table  1  Comparison  of  numerical  and  experimental 
measurements  of  iw  ID  for  ReD  -  9500. 


Time 

Numerical 

Result 

Experimental 

Result 

0.50 

0.025 

0.028 

0.75 

0.061 

0.069 

1.00 

0.098 

0.097 

1.25 

0.135 

0.167 

1.50 

0.232 

0.250 

1.75 

0.293 

0.389 

2.00 

0.440 

0.528 

Table  2  Comparison  of  numerical  and  experimental 
measurements  of  a^J  D  for  ReD  =  9500. 


Time 

Numerical 

Result 

Experimental 

Result 

0.75 

-0.082 

-0.094 

1.00 

-0.061 

-0.044 

1.25 

0.013 

0.013 

1.50 

0.071 

0.069 

1.75 

0.075 

0.113 

Table  3  Comparison  of  numerical  and  experimental 
measurements  of  /  2D  for  ReD  =  9500 . 


Time 

Numerical 

Result 

Experimental 

Result 

0.75 

0.348 

0.350 

LOO 

0.427 

0.363 

1.25 

0.348 

0.300 

1.50 

0.285 

0.289 

1.75 

0.329 

0.325 

of  the  medium  grid;  a  similar  25%  increase  of  the 
medium  grid  is  used  to  obtain  a  fine  (881  x  113)  grid 
The  results  obtained  on  each  of  the  aforementioned 
grids  are  quantitatively  examined  in  terms  of  the 
nondimensionalized  wake  length  streamwise  and 
normal  directions  lw!D.  Additionally,  some 
discussion  is  included  regarding  the  qualitative 
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differences  observed  in  the  flow  structures  obtained 
using  the  different  grids. 

Before  presenting  the  simulation  results,  it  is 
necessary  to  give  some  information  pertaining  to  the 
standards  used  to  assert  grid  independence  of  the 
solutions,  since  a  vast  number  of  methods  have  been 
developed  by  researchers  in  the  reporting  of  grid 
refinement  studies.  In  the  current  study,  the  Grid 
Convergence  Index  (GCI)  of  Roache  (1993)  is  used  in 
analyzing  the  simulation  results.  This  technique  was 
developed  to  provide  a  uniform  method  of  reporting  of 
grid  refinement  studies,  thus  alleviating  the  difficulties 
associated  with  the  varied,  and  sometimes  confusing, 
methods  commonly  found  in  the  literature.  The  use  of 
the  GCI  accomplishes  this  by  relating  the  degree  of 
error  realized  in  considering  a  p-th  order  accurate 
solution  method  and  two  different  grids  of  grid-size 
ratio  r  to  the  degree  of  error  which  would  have  been 
obtained  via  a  grid  refinement  study  of  the  same 
problem  with  the  same  base  grid  using  p-  2  and 
r  =  2 .  Details  of  the  method  used  to  compute  the  GCI 
for  the  fine  and  coarse  grids  is  given  in  Roache  (1993). 


Table  4  Comparison  of  values  of  lw  /  D  for  base  grid, 
coarse  grid  and  experiment. 


Time 

Base  Grid 
Result 

Coarse 
Grid  Result 

Experimental 

Result 

0.50 

0.025 

0.015 

0.028 

0.75 

0.061 

0.049 

0.069 

1.00 

0.098 

0.095 

0.097 

1.25 

0.135 

0.129 

0.167 

1.50 

0.232 

0.186 

0.250 

1.75 

0.293 

0.287 

0.389 

2.00 

0.440 

0.438 

0.528 

Table  5  Comparison  of  values  of  lw  /  D  for  base  gird, 
fine  grid  and  experiment. 


Time 

Base  Grid 
Result 

Fine  Grid 
Result 

Experimental 

Result 

0.50 

0.025 

0.026 

0.028  | 

0.75 

0.061 

0.064 

0.069 

1.00 

0.098 

1.000 

0.097 

1.25 

0.135 

0.136 

0.167 

1.50 

0.232 

0.232 

0.250 

1.75 

0.293 

0.293 

0.389 

2.00 

0.440 

0.445 

0.528 

Examination  of  the  vorticity  and  stream 
function  results  obtained  for  the  coarse  grid  clearly 


showed  the  inability  of  this  grid  to  adequately  resolve 
the  flow  characteristics.  This  is  reflected  in  the 
quantitative  comparison  of  the  measured  values  of 
lw ID  obtained  for  this  case  with  the  base  grid  and 
experimental  results,  given  in  Table  4.  For  the  fine 
grid,  the  flow  structures  were  found  to  appear  very 
similar  to  those  obtained  previously  using  the  medium 
grid,  indicating  that  the  medium  grid  yields  sufficient 
resolution  of  the  flow  structures.  A  comparison  of 
values  of  I D  obtained  for  this  case  with  the  both 
the  medium  grid  and  experimental  results  is  shown  in 
Table  5,  and  a  veiy  close  agreement  is  shown  in  the 
results. 

To  quantitatively  compare  the  grid 
independence  study  results  thus  obtained,  the  Grid 
Convergence  Index  ( GCI)  is  computed  for  both  the 
coarse  and  fine  grids.  These  results  are  shown  in 
Table  6.  It  is  thus  shown  that  the  coarse  grid  yields 
insufficient  resolution  of  the  flow  structures  and  hence 
does  not  yield  a  grid  independent  result.  However,  the 
fine  grid  GCI  shows  a  great  improvement  and  thus 
demonstrates  the  grid  independence  of  the  results 
obtained  using  the  base  and  fine  grids.  Note  that  both 
the  fine  and  coarse  grid  values  of  the  GCI  generally 
show  improvement  at  later  times.  This  may  be 
attributed  to  the  fact  that  the  initial  condition  is  not 
grid  independent,  since  the  calculation  of  body  surface 
vorticity  values  at  time  t  =  0+  depends  on  the  distance 
of  the  first  normal  grid  point  from  the  body  surface. 

Table  6  Comparison  of  values  of  lw  /  D  for  base  grid, 

fine  grid  and  experiment. 


Time 

Coarse  Grid 
GCI 

Fine  Grid 
GCI 

0.50 

2.213 

0.213 

0.75 

1.002 

0.267 

1.00 

0.202 

0.107 

1.25 

0.180 

0.040 

1.50 

0.741 

0.000 

1.75 

0.021 

0.000 

2.00 

0.075 

0.061 

6.  Conclusions  and  Recommendations 

The  self-consistent  formulation,  using 
circulation  as  a  means  to  arrive  at  the  correct  far-field 
boundary  condition  for  the  stream  function,  is 
presented.  The  unsteady  Navier-Stokes  analysis 
presented  for  the  study  of  incompressible  viscous  flows 
over  Joukowski  airfoils  has  been  validated  utilizing  the 
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circular  cylinder  geometry,  due  to  the  wealth  of 
experimental  and  numerical  data  which  exists  for  this 
flow  configuration.  Therefore,  flow  past  a  circular 
cylinder  was  simulated  for  flow  configurations  for 
which  Bouard  and  Coutanceau  (1980)  have  given 
satisfactory  experimental  data.  A  high  degree  of 
agreement  is  found  between  the  numerical  simulation 
results  and  the  experimental  results  for  both  the  simple 
test  case  of  ReD  -  200  and  also  for  the  more  complex 

flow  case  of  ReD  =  9500 .  Thus,  the  self-consistent 

formulation  utilized  here  permits  an  alternative 
implementation  of  the  far-field  boundary  condition  for 
the  disturbance  stream  function,  while  still  retaining 
the  original  surface  boundary  condition  for  the 
vorticity.  Indeed,  the  results  obtained  for  early  time 
integration  show  a  remarkably  high  degree  of 
conformity  with  the  experimental  data  in  the 
resolution  of  the  secondaiy-eddy  phenomenon  for  this 
case.  Following  the  method  of  Yang  (1992),  the  use  of 
upwind-differencing  has  been  successful  in  averting 
the  difficulties  associated  with  the  stability  restrictions 
of  the  central  difference  scheme.  No  spurious 
oscillations  in  the  vorticity  field  were  observed  when 
necessary  grid  sizes  were  used,  even  for  the  higher-Re 
flow  case  of  ReD  =  9500 . 

The  new  formulation  was  also  implemented  in 
an  unsteady  analysis  for  flow  past  a  maneuvering  body 
by  Osswald  at  CFDRL.  His  observations  were  that, 
although  the  force  field  results  are  considerably 
improved  over  the  existing  (a, if/)  formulation,  the 
numerical  results  nevertheless  do  not  conform  to  the 
experimental  results  in  the  highly  nonlinear  unsteady 
flow  regime.  Further  examination  of  the  current 
formulation  suggests  that  the  vorticity  boundary 
condition  at  the  surface  may  need  re-examination.  In 
the  meantime,  until  further  improvements  in  the 
surface  boundary  condition  are  obtained,  the  successful 
validation  thus  realized  for  the  current  formulation 
makes  it  desirable  to  extend  the  present  work  to  study 
more  complex  flow  cases,  including  flows  over 
Joukowski  airfoils  at  various  angles  of  attack. 
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Figure  3.  Schematic  representation  of  measured  wake  properties. 
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Abstract 

A  2-D  unsteady  analysis  is  developed  using  the  Navier- 
Stokes  equations  in  generalized  coordinates.  The  use  of 
generalized  curvilinear  coordinates  permits  flexibility  in 
geometry.  A  spectral  method  is  employed  to  achieve  the 
desired  higher-order  accuracy.  Further,  the  utilization 
of  the  multidomain  technique  facilitates  the  resolution  of 
the  disparate  length  scales  associated  with  Navier-Stokes 
solutions.  The  analysis  is  validated  for  a  model  geome¬ 
try  and  base  flow  results  are  obtained  for  the  parabolic 
cylinder  geometry.  The  ultimate  goal  of  this  work  is  to 
investigate  the  leading-edge  receptivity  to  free-stream 
acoustic  disturbances  on  a  parabolic  cylinder. 


1  Introduction 

Transition  to  turbulence  plays  a  critical  role  in  many  en¬ 
gineering  problems  ranging  from  heat  transfer  to  drag  re¬ 
duction  on  such  applications  as  the  National  Aerospace 
Plane  and  the  High-Speed  Civil  Transport  (HSCT)  ve¬ 
hicle.  Drag-reduction  technologies  are  especially  impor¬ 
tant,  considering  the  saving  that  can  be  obtained  for  the 
commercial  airline  industry  by  achieving  laminar  flow 
on  some  or  all  exterior  surfaces  through  Laminar  Flow 
Control  (LFC)  [1].  LFC  seeks  to  delay/eliminate  tran¬ 
sition  by  attenuating  the  growth  of  instability  modes. 
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Additionally,  in  conjunction  with  wave-  and  vortex-drag 
reduction  techniques,  LFC  can  radically  alter  both  the 
economic  and  ecological  viability  issues  associated  with 
the  HSCT  [2]. 

Receptivity,  the  earliest  stage  of  transition  to  turbu¬ 
lence,  is  the  process  by  which  free-stream  disturbances 
get  internalized  inside  of  a  boundary  layer  in  the  form  of 
instability  waves,  e.g.,  Tollmien-Schlichting  waves.  Re¬ 
ceptivity  complements  linear  stability  theory  in  that  it 
provides  what  has  been  lacking  from  the  theory,  mainly 
the  initial  amplitude  of  the  instability  wave.  The  initial 
amplitude  would  make  the  results  obtained  from  empir¬ 
ical  methods  such  as  the  eN  method  of  Smith  and  Gam- 
beroni  [3]  and  Van  Ingen  [4]  more  reliable.  Currently,  re¬ 
ceptivity  is  the  least  understood  phase  of  the  transition 
process.  Of  primary  concern  is  how  the  long- wavelength 
disturbance  wave  is  able  to  transfer  energy  to  the  much 
shorter-wavelength  instability  wave.  The  high- Reynolds 
number  asymptotic  studies  of  Goldstein  [5,6]  and  Ker- 
schen  [7]  indicate  that  rapid  adjustments  in  the  stream- 
wise  direction  of  the  mean  flow  provide  a  wavelength 
conversion  mechanism.  A  region  where  such  stream  wise 
adjustments  occur  is  the  leading-edge  region,  where  the 
boundary  layer  grows  rapidly. 

Theory  and  experiments  have  revealed  three  possible 
mechanisms  for  the  problem  of  leading-edge  receptivity 
in  the  absence  of  roughness.  In  the  first,  Goldstein's 
leading-edge  receptivity  mechanism  [5],  the  wavelength 
is  shortened  by  the  nonparallel  effects  of  the  rapidly 
growing  boundary  layer.  Goldstein  [5]  developed  this 
theory  by  examining  the  receptivity  of  an  infinitely-thin 
flat  plate  to  a  plane  acoustic  wave.  Heinrich  et  al.  [8] 
later  refined  this  analysis  to  include  oblique  acoustic 
waves  and  convected  gusts.  It  was  found  that  oblique 
acoustic  waves  provide  a  much  higher  receptivity  level 
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than  the  plane  acoustic  wave,  due  to  an  unsteady  flow 
around  the  leading-edge.  In  both  theories,  substantial 
decay  of  the  instability  waves  occurs  prior  to  their  reach¬ 
ing  the  neutral  stability  point  where  they  could  begin  to 
grow.  The  infinitely-thin  flat  plate  might  be  the  only 
model  suitable  for  isolating  this  non-parallel  receptivity 
mechanism,  thus  possibly  excluding  experimental  inves¬ 
tigation  of  this  receptivity  mechanism.  The  numerical 
studies  of  Murdock  [9]  and  Gatski  and  Grosch  [10]  on 
the  infinitely-thin  flat  plate  have  been  inconclusive  due 
to  the  exclusion  of  the  leading-edge  region  (Murdock) 
and  the  receptivity  results  being  preliminary  in  nature 
(Gatski),  i.e.,  early-time  results.  Irregardless,  this  non¬ 
parallel  boundary-layer  receptivity  mechanism  can  be 
considered  of  higher  order  as  compared  to  the  other  two 
mechanisms  described  in  the  following  paragraphs. 

The  second  receptivity  mechanism,  also  formulated  by 
Goldstein  [6],  corresponds  to  the  scattering  of  the  dis¬ 
turbance’s  long  wavelength  on  some  short-scale  stream- 
wise  variation.  With  respect  to  the  leading-edge  recep¬ 
tivity  problem,  the  surface  variation  could  occur  in  the 
form  of  a  surface  curvature  discontinuity  somewhere  in 
the  leading-edge  geometry.  Strong  evidence  of  this  re¬ 
ceptivity  mechanism  is  provided  by  the  experiments  of 
Wlezien  et  al.  [11-13]  and  the  computations  of  Reed  et 
al.  [14, 15]  on  a  semi-infinite  flat  plate  with  an  elliptical 
leading  edge. 

The  third  and  final  leading-edge  receptivity  mecha¬ 
nism  present  is  given  by  Nishioka  and  Morkovin  [16] 
and  relies  on  the  pressure  gradient  field  to  provide  the 
needed  wavelength  conversion  mechanism.  In  essence, 
this  mechanism  operates  in  a  manner  similar  to  that 
of  the  surface- variation  receptivity  mechanism,  with  the 
important  distinction  that  the  short-scale  streamwise 
variations  are  provided  by  the  pressure  field  rather  than 
by  the  surface.  This  receptivity  mechanism  should  be 
active  to  some  degree  for  all  blunt-body  geometries.  Ev¬ 
idence  of  this  mechanism  in  leading-edge  geometries  is 
provided  by  the  numerical  results  of  Reed  et  al.  [14, 15] 
and  Buter  et  al.  [17, 18]  for  vortical  disturbances. 

Lastly,  it  should  be  noted  that  the  inherent  favor¬ 
able  or  adverse  pressure  gradients  associated  with  the 
leading-edge  receptivity  problem  also  lead  to  a  retarda¬ 
tion  or  acceleration,  respectively,  of  the  entire  transition 
process  by  shifting  the  neutral  stability  curve  towards 
or  away  from  the  leading  edge,  respectively. 

In  the  present  study,  the  analysis  and  numerical  proce¬ 
dure  are  developed  for  studying  the  leading-edge  recep¬ 
tivity  on  a  parabolic  cylinder.  The  analysis  is  validated 
using  the  shear-driven  cavity  geometry.  Subsequently,  it 
is  extended  to  the  parabolic  cylinder,  and  steady-state 


base  flows  are  obtained  and  compared  with  the  work  of 
Davis  [19].  In  addition,  the  resolution  requirements  for 
the  receptivity  problem  are  examined. 


2  Mathematical  Formulation 


The  flow  past  a  semi-infinite  parabolic  cylinder  is  simu¬ 
lated  by  employing  the  2-D,  unsteady  vorticity /stream 
function  (u>,V0  form  of  the  Navier-Stokes  equations  in 
generalized  curvilinear  coordinates: 

Vorticity  Transport  Equation: 

-du  3  (  drl>\  d  (  drj>\ 

^  dt  +  ae  vae)  de  v"aev  “ 
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Stream  Function  Equation: 
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where  xp  is  the  total  stream  function  and  ipv  is  the  vis¬ 
cous  deviation  from  the  inviscid  stream  function, 
i.e.,  xp  =  rpinv  -f  xpv .  The  viscous  deviational  stream 
function  has  been  introduced  to  alleviate  the  unbound¬ 
edness  of  ip  by  analytically  containing  it  in  xpinv .  In  ad¬ 
dition,  (f1,^2)  are  the  generalized  orthogonal  curvilin¬ 
ear  coordinates,  with  metrics,  y/g ,  </n,  <722,  determined 
by  the  transformation  between  the  physical  coordinates 
and  the  generalized  coordinates.  The  equations  have 
been  nondimensionalized  using  the  kinematic  viscosity, 
j/,  the  free-stream  velocity,  [/«>,  and  the  frequency  of 
the  free-stream  disturbance,  f),  to  produce  the  Reynolds 

vt 


number  Re  = 


This  nondimensionalization  is  con¬ 


sistent  with  that  used  by  Goldstein  [5,6]. 


2.1  Boundary  Conditions 

The  appropriate  boundary  conditions  consist  of  the  no¬ 
slip  and  no-penetration  conditions  on  the  body  surface, 
free-stream  conditions  in  the  far-field,  symmetry  condi¬ 
tions  along  the  stagnation  line,  and  outflow  conditions 
corresponding  to  the  Blasius  boundary-layer  solution 
at  downstream  infinity.  Mathematically,  the  boundary 
conditions  can  be  stated  as  follows 

w  =  xpv  =  0  along  symmetry  line, 
dxpv 

w  =  -r-r-  =0  at  far-field  infinity, 

d £2 
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on  the  body  surface . 


(3) 


2,2  Solution  Procedure 

The  overall  solution  procedure  is  subdivided  into  two 
phases:  the  determination  of  the  base  flow  and  the  sim¬ 
ulation  of  its  receptivity.  For  the  base  flow,  the  work 
of  Davis  [19]  has  served  as  a  guide.  The  domain  of  in¬ 
terest  was  reduced  by  exploiting  the  symmetry  of  the 
flow  about  the  stagnation  line.  Davis’  nondimesional- 
ization  was  also  adopted  for  this  phase  of  the  analy¬ 
sis.  The  Reynolds  number,  Rer  =  ,  is  based  on 

the  nose  radius  rn,  and  the  independent  spatial  coor¬ 
dinates  are  scaled  according  to  Rer.  In  addition,  the 
similarity  variables  introduced  by  Davis  were  utilized  to 
aid  in  the  application  of  the  outflow  boundary  condi¬ 
tions.  Specifically,  the  similarity  variables  were  used  in 
the  farthest  downstream  domain  to  allow  the  Blasius  so¬ 
lution  to  be  applied  accurately.  The  similarity  variables 
were  not  used  in  the  domain  around  the  leading  edge 
as  the  similarity-variable  form  of  the  governing  equa¬ 
tions  contain  terms  which  are  singular  at  the  stagnation 
line.  Spectral  methods  are  notorious  for  being  sensitive 
to  singularities  and  this  singularity  was  no  exception. 
Since  a  steady-state  base  solution  is  desired,  a  fictitious 
time-derivative  was  added  to  Eq.  (2)  to  accelerate  the 
convergence  to  the  steady-state  solution. 

For  the  receptivity  simulation,  the  base  flow  is  re- 
nondimensionalized  to  be  consistent  with  Eqs.  (1-2). 
The  plane  acoustic  wave  is  simulated  by  pulsating  the 
free-stream  flow,  i.e., 

Uoo(t)  =  1  +  csint ,  (4) 

where  e  is  a  small  quantity. 


3  Numerical  Method 

The  algebraic  equations  are  obtained  by  discretizing 
spectrally  in  space  and  using  a  finite-difference  approx¬ 
imation  in  time  [20].  A  spectral  collocation  method  has 
been  chosen  over  other  numerical  differencing  schemes 
due  to  its  high  accuracy  and  negligible  phase  errors. 
Both  of  these  characteristics  will  be  needed  to  simulate 
the  small  amplitude  waves.  The  dependent  variables  are 
expanded  using  Chebyshev  polynomials  in  both  space 
dimensions,  and  the  governing  equations  are  enforced 
at  the  Gauss-Lobatto  collocation  points.  All  diffu¬ 
sive  terms  are  treated  implicitly  using  a  Crank-Nicolson 
scheme,  while  the  nonlinear  convective  terms  are  ad¬ 
vanced  using  a  second-order  Adams-Bashforth  scheme. 


The  resulting  set  of  semi-implicit  algebraic  equa¬ 
tions  is  solved  iteratively  at  each  time  step  us¬ 
ing  a  Preconditioned  Truncated  Conjugate  Residual 
(PTCR)  method  [21],  Finite-difference  and  finite- 
element  preconditioners  are  employed  for  the  time- 
dependent  and  time-independent  equations,  respec¬ 
tively.  The  finite-difference  preconditioners  are  in¬ 
verted  by  an  Alternating-Direction  Implicit  method, 
while  the  finite-element  preconditioner  uses  an  efficient 
block  gaussian  elimination  scheme. 

Since  spectral  methods  are  extremely  sensitive  to  the 
boundary  conditions  used,  care  has  been  exercized  in  the 
implementation  of  the  boundary  conditions  so  as  to  not 
degrade  the  accuracy  of  the  method.  Particular  atten¬ 
tion  has  been  paid  to  the  no-slip  condition  which  was 
implemented  via  an  integral  constraint  condition  [22]. 

In  addition,  the  spectral  method  is  supplemented  by 
a  patched  multidomain  technique  which  ensures  C1- 
continuity  on  the  interface.  Unfortunately,  the  bet¬ 
ter  global-integral-flux  interface  conditions  of  Macaraeg 
et  al.  [23]  could  not  be  used  due  to  the  similarity- 
variable  form  of  the  equations  employed  in  the  last 
subdomain.  The  interface  conditions  and  the  surface 
boundary  conditions  have  been  implemented  using  the 
influence-coefficient-matrix  technique  [23].  Three  sub- 
domains  were  employed  in  the  streamwise  direction.  The 
first  two  allow  adequate  resolution  of  the  region  of  in¬ 
terest,  while  the  remaining  domain  served  as  a  buffer 
domain  for  the  receptivity  calculations. 

The  grid  is  generated  using  parabolic  coordinates  in 
conjunction  with  1-D  clustering  transformations.  The 
clustering  transformations,  along  with  the  multiple  do¬ 
mains,  allow  the  disparate  length  scales  of  the  Navier- 
Stokes  solutions  to  be  resolved.  For  the  parabolic  cylin¬ 
der,  the  use  of  multiple  domains  in  the  streamwise  di¬ 
rection  also  facilitates  the  implementation  of  the  buffer 
domain  technique  [24]  in  which  the  small  disturbances 
are  passed  out  of  the  buffer  domain  without  reflection. 
This  takes  the  place  of  the  outflow  boundary  condition 
for  the  case  in  which  the  acoustic  disturbance  has  been 
introduced.  The  far-field  boundary  is  located  at  a  large, 
but  finite  normal  distance  from  the  body  surface.  Care 
has  been  exercized  to  ensure  that  this  location  is  suffi¬ 
ciently  far  away  so  as  to  not  influence  the  solution.  In 
the  streamwise  direction,  the  grid  generation  has  em¬ 
ployed  the  logarithmic  stretching,  used  by  Davis  [19], 
in  the  last  domain  to  stretch  to  downstream  infinity. 
This  allows  for  the  accurate  implementation  of  the  Bla¬ 
sius  solution  as  the  outflow  boundary  condition  for  the 
base  flow.  The  spectral  method  experienced  no  adverse 
effects  since  the  similarity  variables  were  used  in  this 
last  domain.  In  addition,  the  parameters  of  the  normal 
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stretching  transformation  were  set  such  that  half  of  the 
total  number  of  grid  points  in  the  normal  direction  lay 
within  the  boundary  layer. 


4  Results 

The  overall  analysis  and  the  numerical  technique  devel¬ 
oped  were  first  validated  using  the  results  obtained  for  a 
shear-driven  square  cavity  for  various  Reynolds  numbers 
and  compared  successfully  with  the  available  benchmark 
solutions  of  Ghia  et  al.  [25]  and  Nishida  et  al.  [26].  Ta¬ 
ble  1  shows  a  comparison  of  the  primary  vortex  prop¬ 
erties  for  Re  =  100  as  obtained  by  the  present  method 
for  various  discretizations  with  both  Ghia  et  al.’s  [25] 
2nd-order  accurate  results  and  Nishida  et  al.’s  [26]  10th- 
order  accurate  results  on  129  x  129  grids.  The  results 
for  the  present  method  have  been  obtained  on  three 
grid  discretizations  using  both  Chebyshev  and  Legendre 
polynomials  as  expansion  bases.  Excellent  agreement 
is  achieved  despite  the  coarse  discretization  used  in  the 
present  method.  Figure  1  depicts  some  typical  plots 
of  the  grid,  stream  function,  vorticity,  and  vorticity  for 
Re  =  100. 

Subsequently,  results  are  obtained  for  flow  past  a 
parabolic  cylinder. 

4.1  Base  Flow 

Computations  on  a  single  domain  were  performed  to  as¬ 
certain  the  effects  of  the  placement  of  the  normal  far- 
field  boundary  on  the  solution,  as  well  as  determine  the 
grid  size  necessary  for  accurate  resolution  of  the  prob¬ 
lem.  It  was  found  that  the  placement  of  the  far-field 
boundary  had  little  effect  on  the  solution  for  a  range  of 
values  of  the  distance  of  this  boundary  from  the  body 
surface.  In  addition,  essentially  similar  results  are  ob¬ 
tained  with  discretizations  as  small  as  17  x  25  points. 

Figure  2  depicts  the  wall  vorticity  function  for  various 
Reynolds  numbers  for  the  three-domain  configuration. 
The  results  were  obtained  using  a  total  of  53  points  in 
the  streamwise  direction  for  all  three  domains  and  33 
points  in  the  normal  direction.  There  is  no  discernable 
indication  in  the  solution  as  to  the  presence  of  the  inter¬ 
faces  between  the  domains.  The  results  shown  in  Fig.  2 
conform  to  those  obtained  by  Davis  [19]  which  are  re¬ 
produced  here  in  Fig.  3  (The  reader  is  asked  to  ignore 
the  symbols  which  correspond  to  the  parabolized  Navier- 
Stokes  solutions).  Davis  used  a  finite-difference  method 
with  a  40  x  125  grid.  In  each  of  the  present  calcula¬ 
tions,  the  interface  between  the  last  two  domains  lies  in 


the  region  where  the  wall  vorticity  function  has  essen¬ 
tially  attained  the  Blasius  value  (See  Fig.  2).  This  is  so 
that,  in  the  receptivity  calculations,  the  instability  wave 
will  have  a  chance  to  grow,  i.e.,  the  flow  will  pass  branch 
I  of  the  linear  stability  curve  prior  to  entering  the  buffer 
domain. 

To  show  that  the  flow  variables  are  continuous  across 
the  multidomain  interfaces,  Fig.  4  depicts  the  compu¬ 
tational  plane  for  the  three-domain  configuration,  along 
with  the  contours  of  the  similarity  viscous-deviational 
stream  function,  similarity  vorticity  and  physical  vortic¬ 
ity  for  Rer  =  100.  As  can  be  clearly  seen,  there  is  no 
discontinuity  in  the  flow  variables  across  the  interfaces, 
despite  the  discontinuity  in  grid  metrics.  Figure  5  shows 
some  physical-plane  results  for  Rer  =  100;  the  contours 
here  are  of  the  total  stream  function  and  vorticity.  The 
extent  of  the  physical  domain  shown  in  this  figure  is 
contained  solely  in  the  first  subdomain. 


4.2  Receptivity  Requirements 

Before  embarking  on  the  simulation  of  the  receptivity, 
it  is  prudent  to  examine  what  has  been  done  previously 
and  what  the  grid  requirements  will  be  for  accurately 
resolving  the  Tollmien-Schlichting  wavelength. 

Murdock  [27]  first  examined  this  geometry’s  response 
to  free-stream  acoustic  disturbances  using  the  Parabo¬ 
lized  Navier-Stokes  equations.  His  results  indicated  a 
declining  receptivity  level  with  increasing  nose  radius. 
Unfortunately,  he  had  difficulties  with  his  downstream 
boundary  condition  [28].  More  recently,  Kerschen  and 
co-workers  [28,29]  and  Gatski  and  Kerschen  [30]  have 
collaborated  to  study  this  geometry  using  theoretical 
and  computational  tools,  respectively.  The  theoretical 
work  of  Hammer  ton  et  al.  [29],  based  on  a  matched- 
asymptotic  expansion  analysis,  has  also  indicated  that 
the  effect  of  nose  bluntness  is  to  decrease  the  receptivity. 
They  introduced  a  Strouhal  number,  S  = 
which  is  a  measure  of  the  nose  bluntness,  'fke  recep¬ 
tivity  was  found  to  decrease  by  an  order  of  magnitude 
when  S  increased  to  a  value  of  0.3.  Gatski  et  al.  [30] 
presented  some  preliminary  results  for  one  receptivity 
case.  His  disturbance  quantities  compare  qualitatively 
with  the  Stokes-shear-wave  response  to  acoustic  forcing. 

As  a  further  validation  of  the  present  method,  it 
is  planned  to  simulate  the  same  receptivity  case  that 
Gatski  et  al.  [30]  examined.  The  parameters  for  this  case 
are  as  follows:  Re  —  104  and  Rer  =  10.  The  nose  ra¬ 
dius  for  this  case  is  smaller  than  the  Tollmien-Schlichting 
wavelength  and  the  favourable  pressure  gradient  is  not 
as  strong  as  for  larger  nose  radii.  Therefore,  based  on 
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the  theory  [29],  the  receptivity  can  be  expected  to  be 
nearly  optimal  for  the  configuration.  The  Tollmien- 
Schlichting  wavelength  in  the  current  nondimensional- 
ization  will  be  approximately  2  [30].  Spectral  methods 
need  at  least  three  points  per  wavelength,  thus  requir¬ 
ing  approximately  100  points  or  more  in  the  streamwise 
direction  distributed  in  the  first  two  subdomains. 

Of  the  three  basic  receptivity  mechanisms  discussed 
in  the  introduction,  only  Goldstein’s  [5]  leading-edge 
and  Nishioka  et  al.’s  [16]  receptivity  mechanisms  can 
be  considered  active  for  the  parabolic  cylinder  recep¬ 
tivity  problem  since  the  parabolic  cylinder  geometry  is 
free  of  surface  variation  discontinuities.  In  addition,  the 
nose  bluntness  produces  a  favourable  pressure  gradient 
over  a  portion  of  the  nose  region,  further  retarding  the 
receptivity  process.  Hence,  any  receptivity  can  be  ex¬ 
pected  to  be  small  compared  to  that  of  other  receptivity 
mechanisms.  This  places  a  severe  strain  on  any  numer¬ 
ical  method  used  to  simulate  this  receptivity  problem. 
These  results  are  presently  being  generated. 


5  Conclusion 

An  analysis  and  a  spectral  multidomain  method  are  de¬ 
veloped  for  examining  receptivity  problems.  The  anal¬ 
ysis  and  method  are  validated  using  a  model  geometry, 
and  thereafter,  extended  to  the  parabolic  cylinder  geom¬ 
etry.  Steady-state  base  flows  are  obtained  and  conform 
to  previous  numerical  results.  Lastly,  the  requirements 
for  simulating  receptivity  to  plane  acoustic  waves  are 
examined. 
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Figure  3:  Solutions  for  Wall  Vorticity  Function  by  Davis  [19]. 
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Figure  4:  Parabolic  Cylinder  Results  Plotted  in  Computational  Plane:  Rer  =  100. 
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Figure  5:  Physical  Plane  Results  for  Rer  =  100. 
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ABSTRACT 

Leading-edge  receptivity  to  tree-stream  acoustic 
disturbances  is  studied.  A  semi-infinite  geometry,  i.e.,  the 
parabolic  cylinder,  is  employed  to  isolate  the  effects  of  nose- 
bluntness  on  the  receptivity  process.  The  receptivity  process  is 
simulated  using  the  2-D  unsteady  Navier-Stokes  equations  in 
generalized  curvilinear  coordinates.  In  addition  to  the 
generalized  coordinates,  a  multidomain  technique  is  utilized  to 
help  resolve  the  disparate  length  scales  of  the  problem.  A 
spectral  collocation  method  has  been  selected  for  the  spatial 
discretization  since  the  receptivity  mechanisms  active  in  this 
geometry  are  weak  and  the  resulting  instability  wave  can  be 
expected  to  undergo  substantial  decay  prior  to  reaching  the 
neutral  stability  location.  The  overall  method  is  validated  by 
comparing  the  steady-state  base  flow  results  with  previously 
published  numerical  results.  Results  for  Reynolds  numbers  of 
10,  100,  and  1000  are  obtained  where  the  Reynolds  number  is 
based  on  the  nose  radius  of  the  parabola.  This  allows 
assessment  of  the  effects  of  nose-bluntness  on  the  receptivity 
mechanism.  Forcing  is  introduced  by  pulsating  the  free- stream 
velocity  to  simulate  an  acoustic  wave.  Both  plane  and  oblique 
acoustic  waves  can  be  simulated  in  this  manner.  The  forcing 
frequency  has  been  selected  such  that  the  Reynolds  number 
based  on  the  forcing  frequency  is  10,000. 

1.  INTRODUCTION 

Transition  to  turbulence  continues  to  be  of  critical 
importance  to  a  variety'  of  current  technologies  such  as  drag 
reduction  and  heat  transfer.  Several  conceptual  models  of  the 
many  routes  to  transition  exist  [1],  but  the  most  prevalent  model 
is  one  in  which  the  transition  process  is  divided  into  three 
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phases:  receptivity,  linear  stability',  and  nonlinear  breakdown. 
The  linear  stability  phase  has  been  well  established  for  many 
years  and  has  provided  the  current  tools  used  in  transition 

prediction,  namely,  the  e  method  [2].  Similarly,  the  nonlinear 
breakdown  phase,  while  not  completely  understood,  is 
represented  fairly  well  in  the  early  stages  by  secondary 
instability  theory  [3].  The  first  phase,  receptivity,  is  the  process 
by  which  small  disturbances  in  the  free  stream  enter  the 
boundary  layer  and  excite  the  instability'  modes  of  the  boundary 
layer.  This  phase  of  transition  has  been  the  least  understood; 
however,  in  recent  years,  significant  progress  has  been  made 
towards  its  theoretical  understanding.  What  has  been  lacking  is 
an  understanding  of  how  free  stream  disturbances,  typically  of 
much  longer  wavelength  than  instability'  waves,  are  able  to 
transfer  energy  into  instability  wave  in  the  presence  of  a 
wavelength  mismatch.  Theories  by  Goldstein  [4]  and  Kerschen 
[5]  indicate  that  rapid  adjustments  in  the  stream  wise  direction 
provide  the  necessary  wavelength  conversion  mechanism 
resulting  in  the  required  wavelength  match-up.  Regions  where 
such  rapid  adjustments  occur  include  the  leading-edge  region 
where  the  boundary'  layer  is  growing  rapidly. 

In  the  present  study,  leading-edge  receptivity  to  a  free- 
stream  acoustic  wave  is  examined  for  a  blunt  leading-edge 
geometry  consisting  of  a  parabolic  cylinder.  The  primary  goals 
of  this  work  are  to  proride  verification  of  some  recent  theoretical 
results  obtained  by  Hammerton  et  al.  [6]  for  the  effects  of  nose 
bluntness  on  receptivity  and  to  provide  a  numerical  tool  for 
studying  receptivity  for  flows  in  non-simple  geometries. 

2.  MATHEMATICAL  AND  NUMERICAL  METHODS 

In  the  present  analysis,  the  2-D  unsteady  Navier- 
Stokes  equations  are  utilized  in  the  vorticity/stream  function 
formulation.  In  addition,  the  equations  have  been  developed  in 
generalized  curvilinear  coordinates  to  better  represent  arbitrary 
geometries  and  to  achieve  improved  resolution  of  the  disparate 
length  scales  of  the  flow  physics  Following  the  theory  on 
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receptivity,  the  equations  have  been  nondimensionalized  using 
the  free-stream  velocity,  the  free-stream  kinematic  viscosity  and 
the  forcing  frequency  of  the  free-stream  disturbance.  This 
provides  the  appropriate  scaling  of  the  variables  for  the 
receptivity  study.  The  necessary  boundary  conditions  are 
provided  by  the  free-stream  conditions  far  from  the  body,  the  no¬ 
slip  and  no-penetration  conditions  at  the  surface,  and  an 
appropriate  outflow  boundary  condition.  Additionally,  the 
symmetry  of  the  problem  is  exploited  to  reduce  the  physical 
domain  of  interest;  hence,  symmetry  boundary  conditions  are 
also  used.  More  details  about  the  analysis  have  been  described 
by  Blodgett  et  al.  [7]. 

A  spectral  collocation  method  [8]  has  been  selected  for 
the  spatial  discretization  because  of  its  low  phase  error  and 
superior  accuracy  properties.  The  equations  are  advanced  in 
time  using  a  semi-implicit  time  stepping  method.  The  resulting 
set  of  algebraic  equations  are  solved  iteratively  using  a 
preconditioned  conjugate  residual  method  [9]  with  a  finite- 
difference  preconditioner. 

In  addition  to  the  spectral  method,  a  multidomain 
method  has  been  added  to  allow  a  zonal  approach  to  the 
problem.  The  utilization  of  the  multidomain  method  allows 
better  use  of  grid  points,  different  grid  clustering  transformations 
in  different  subdomains,  and  even  the  use  of  different  equations 
in  different  subdomains.  The  subdomains  are  connected  by  a 

patched  grid,  and  C1 -continuity  of  the  dependent  variables  is 
enforced  across  the  interfaces  between  the  subdomains. 

The  numerical  implementation  of  the  boundary 
conditions  is  straightforward,  with  the  exception  of  the  no-slip 
condition.  For  this  boundary  condition,  an  integral  constraint 
condition  [10]  has  been  employed.  In  essence,  this  integral 
constraint  condition  allows  the  surface  vorticity  to  be  determined 
so  as  to  be  consistent  with  the  no-slip  condition.  The  remaining 
boundary  conditions  consist  of  symmetry  conditions  along  the 
stagnation  line,  free-stream  conditions  at  the  outer  normal 
boundary,  and  an  outflow  boundary  condition  downstream.  The 
outflow  boundary  condition  will  be  addressed  further  in  the  next 
paragraph.  Lastly,  since  spectral  methods  are  sensitive  to  any 
errors  in  boundary  conditions,  care  has  been  taken  in  the 
implementation  of  the  boundary  conditions. 

For  the  parabolic  cylinder  geometry,  the  benchmark 
work  of  Davis  [11]  has  been  used  to  guide  the  current  analysis. 
His  asymptotic  analysis  of  the  problem  indicates  that  similarity- 
type  variables  are  the  preferred  dependent  variables  for  this 
semi-infinite  geometry.  Unfortunately,  the  resulting  governing 
equations  are  singular  at  the  stagnation  line,  and  spectral 
methods  are  extremely  sensitive  to  singularities.  However,  the 
use  of  the  similarity  variables  provides  an  accurate  outflow 
boundary  condition,  namely,  the  Blasius  boundary-layer  solution. 
Hence,  a  unique  feature  of  the  current  analysis  is  the  use  of 
different  dependent  variables  in  different  subdomains.  In 
particular,  the  similarity  variables  are  utilized  in  the  farthest 
downstream  subdomain  in  the  stream  wise  direction.  This  allows 


the  asymptotically  correct  outflow  boundary  condition  to  be 
applied.  No  difficulty  is  experienced  in  the  solution  process  in 
this  subdomain  since  the  similarity  variables  are  well  behaved  as 
they  approach  downstream  infinity. 

The  grid  is  generated  by  the  conformal  parabolic 
coordinates  as  described  by  Davis  [11].  For  the  present  analysis, 
the  grid  in  the  physical  plane  is  stretched  to  downstream  infinity 
in  the  streamwise  direction,  but  only  to  a  far  but  finite  location 
in  the  normal  direction.  A  three  subdomain  set-up  is  used  in  the 
streamwise  direction,  with  different  clustering  transformations 
possible  in  each  subdomain.  In  particular,  the  streamwise 
clustering  transformation  used  by  Davis  [11]  is  utilized  in  the 
last  subdomain,  while  algebraic  clustering  transformations  are 
used  in  the  first  two  subdomains.  In  the  normal  direction,  a 
clustering  transformation  similar  to  that  used  by  Davis  [11]  is 
employed,  and  allows  the  stretching  parameters  to  be  chosen 
such  that  half  of  the  grid  points  in  the  normal  direction  lie  within 
the  boundary  layer.  The  location  of  the  outer  boundary  in  the 
normal  direction  has  been  selected  carefully  so  that  its  location 
does  not  influence  the  solution  significantly. 

The  receptivity  simulation  is  done  in  two  parts.  In  the 
first  part,  a  base  flow  is  obtained  in  the  absence  of  any  free- 
stream  disturbances.  A  typical  base  flow  might  consist  of  a 
steady-state  solution  to  the  problem.  Once  this  base  flow  is 
obtained,  the  actual  receptivity  simulation  can  commence  by  the 
introduction  of  the  free-stream  disturbance. 

In  obtaining  the  base  flow,  a  fictitious  time  derivative 
has  been  added  to  the  stream  function  equation,  and  local  time 
stepping  is  employed  to  accelerate  the  convergence  of  the  flow 
field  towards  a  steady  state  base  flow.  Despite  the  convergence 
acceleration  techniques,  the  relaxation  time  to  reach  a  steady 
state  is  prohibitively  large  for  the  current  numerical  method. 
Therefore,  once  the  convergence  to  steady-state  reaches  a  certain 
specified  level,  the  solution  is  stopped,  and  this  result  is  used  as 
the  base  flow. 

For  the  receptivity  simulation,  the  time-accurate  form 
of  the  equations  is  employed.  The  free-stream  acoustic 
disturbance  is  introduced  in  the  far  field  by  pulsating  the  free- 
stream  velocity.  The  introduction  of  the  disturbance  poses  an 
additional  difficulty  for  implementing  the  outflow  boundary 
condition.  To  circumvent  this,  the  buffer  domain  technique  [12] 
is  utilized  in  the  last  domain,  so  as  to  attenuate  any  disturbances 
to  zero  in  the  last  domain  without  causing  any  reflection  of 
waves. 

The  present  numerical  method  has  been  validated  on 
the  shear-driven  cavity  geometry  as  well  as  via  the  base  flow 
results  obtained  [7]. 

3.  BASE  FLOW  RESULTS 

Base  flow  results  are  obtained  for  several  parabolas 
with  different  nose  radii.  The  amount  of  nose  bluntness  of  the 
different  parabolas  is  distinguished  by  the  Reynolds  number 
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teased  on  the  nose  radius,  i.e.,  Rey.  Figure  1  depicts  the 
similarity  vorticity  function  along  the  wall  as  a  function  of  the 

stream  wise  similarity  coordinate,  for  various  Re^ .  These 
jesults  compare  very  well  with  those  obtained  by  Davis  [11] 
using  a  finite-difference  method  with  a  40  X  125  grid,  and 
reproduced  here  in  Figure  2.  The  current  results  were  obtained 
with  a  total  of  53  and  33  points  in  the  stream  wise  and  normal 
directions,  respectively.  The  results  shown  in  Fig.  1  have  not 
been  completely  converged  to  a  steady-state  solution,  but  are 
sufficiently  close. 

Figure  3  shows  the  computational  grid,  the  similarity 
vorticity,  the  similarity  viscous-deviational  stream  function,  and 
the  physical  vorticity  in  the  computational  plane  for 

Re  =100.  It  is  clear  from  this  figure  that  the  flow  variables 

r 

are  continuous  across  the  interfaces  which  are  shown  in  the 
contour  plots  as  vertical  lines. 

4.  RECEPTIVITY  SIMULATION 

Receptivity  results  for  three  different  nose  radii,  e.g.. 

Re  =10,  100  and  1000  are  currently  being  obtained  and 

analyzed.  The  forcing  frequency  has  been  selected  to  yield  a 
Remolds  number  of  10,000,  based  on  the  frequency.  The 

amplitude  of  the  free-stream  acoustic  wave  has  been  set  at  1 0-3 
of  the  free-stream  velocity.  This  should  be  sufficiently  small  to 
avoid  nonlinear  responses.  Gatski  and  Kerschen  [13]  have  also 
considered  these  same  parameters  and  have  estimated  that  the 
resulting  Tollmien-Schlichting  instability  wave  should  have  a 
wavelength  around  2.3  units  in  the  current 
nondimensionalization.  In  order  to  resolve  the  Tollmien- 
Schlichting  wavelength,  roughly  three  grid  points  are  required 
per  wavelength  for  the  current  spectral  method.  To  meet  this 
requirement,  the  stream  wise  grid  distribution  has  been  increased 
to  slightly  more  that  150  grid  points  distributed  to  the  first  two 
subdomains.  In  addition,  the  normal  distribution  has  been 
increased  to  49  points.  Since,  a  spectral  method  is  being  used, 
the  base  flow'  results  from  a  coarser  grid  can  be  interpolated  onto 
the  finer  grid  without  loss  of  accuracy. 

Obtaining  results  for  three  different  nose  radii  should 
allow  the  qualitative  reproduction  of  the  trend  predicted  by 
theory  in  Hammerton  et  al.’s  [6]  work.  They  show  that,  as  the 

o  Re 

Strouhal  number,  S  = - L,  is  increased,  the  receptivity 

Rc 

coefficient  decreases. 
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Figure  1:  Navier-Stokes  Solutions  for  Wall  Vorticity  Function  using  Single  Domain. 
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Figure  2:  Solutions  for  Wall  Vorticity  Function  by  Davis  [11]. 
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NUMERICAL  SIMULATION  OF  NONLINEAR  INSTABILITIES 
IN  COMPLEX  CHANNEL  FLOWS* 
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A  tool  for  the  numerical  simulation  of  transition  to  turbulence  in  3D  complex  channel  flows  has  been 
developed  using  a  single-domain  spectral  collocation  method.  The  present  study  describes  results  which 
validate  the  accuracy  of  the  spectral  method  and  demonstrate  the  ability  of  the  code  to  simulate  temporal 
stability  problems.  Nonlinear  stability  results  for  plane  and  curved  channel  flows  were  obtained.  In 
particular,  a  quasi-saturated  Tollmien-Schlichting  wave  for  plane  channel  flow  with  Reynolds  number 
4000  was  computed  and  compares  favorably  with  available  results.  The  nonlinear  evolution  equations 
for  subcritical  plane  and  curved  channel  flows  were  obtained  and  compared  with  Landau’s  equation 
governing  the  growth  of  weakly  nonlinear  waves. 


Introduction 

The  study  of  the  transition  of  laminar  flow  to  a 
turbulent  state  has,  in  recent  research  efforts, 
focused  on  the  nonlinear  growth  and  interaction  of 
waves  which  originate  as  linear  eigenmodes  of  the 
flow-governing  equations.  Many  experimentally  and 
numerically  observed  routes  to  turbulence  can  be 
understood  theoredcally  in  terms  of  the  nonlinear 
instability  of  waves  which  are  excited  in  the  flow. 
Plane  channel  flow  has  been  studied  extensively1-*  and 
has  shown  many  commonalities  with  plane  boundary- 
layer  flow.  In  plane  channel  flow,  Tollmien- 
Schlichting  waves  interact  with  streamwise  vortices  to 
produce  three-dimensional  hairpin,  or  A  vortex, 
structures  which  are  precursors  to  a  fully-turbulent 
state.  Curved  channel  flow,  that  is,  azimuthal  flow  in 
an  annulus,  has  also  been  examined  and  has 
demonstrated  a  different  instability  mechanism  due  to 


streamwise  curvature9-13.  Streamwise  vortices,  called 
Dean  vortices  for  curved  channel  flow,  grow  to  a 
stable  equilibrium  for  a  range  of  Reynolds  numbers. 
Increasing  the  Reynolds  number  in  time  leads  to 
bifurcations  to  more  complex,  but  stable,  states  until 
a  critical  Reynolds  number  is  reached  for  which  no 
stable  flow  exists,  and  the  flow  becomes  fully 
turbulent 

Studies  of  flows  in  more  complex  types  of 
channels  have  shown  that  the  transition  mechanisms 
observed  in  plane  and  curved  channel  flows  are  also 
dominant  in  other  types  of  flows.  Examples  -of 
complex  channel  flows  include  wavy  converging- 
diverging  channel14,  constricted  channel15,16,  grooved 
channel17,1*,  plane  channel  with  2D  and  3D  surface 
perturbations19,  and  channel  flow  with  compliant 
walls20.  In  complex  channel  flows,  the  presence  of 
multiple  shear  layers  at  different  scales  causes  the 
stability  properties  of  the  fundamental  Tollmien- 
Schlich ting-type  modes  and  streamwise  vortex  modes 
to  be  significantly  altered.  For  example  in  the  2D 
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periodic  converging-diverging  channel,  a  chaotic  state 
is  achieved  through  bifurcation  of  the  Tollraien- 
Schlichting  wave  to  a  quasi-periodic  state  as  Reynolds 
number  is  increased,  following  the  Ruelle-Takens- 
Newhouse  scenario14. 

The  purpose  of  the  current  work  is  to  develop  a 
numerical  simulation  tool  to  study  the  temporal 
evoludon  of  instabilities  in  incompressible  channel 
flows  of  various  types,  in  order  to  examine  the  routes 
to  turbulence,  which  are  dominant  in  many  types  of 
flows  of  interest  to  fluid  dynamics  researchers  and  the 
aerospace  industry.  A  numerical  approach  has  been 
developed  for  treating  3D  flows,  forced  by  a  mean 
streamwise  pressure  gradient,  in  arbitrary  channel 
geometries.  These  flows  are  assumed  to  be  periodic  in 
the  streamwise  and  spanwise  directions;  thus,  the 
temporal  evolution  of  disturbances  present  in  the 
initial  flow  is  studied.  Validation  of  the  numerical 
approach  has  included  error  analysis  of  the 
discretization  and  comparison  with  linear  theory  and 
nonlinear  stability  results  from  other  authors.  This 
paper  discusses  the  first  stage  of  the  work,  namely,  the 
validation  of  the  numerical  approach  and  a  study  of 
the  nonlinear  evolution  of  Tollmien-Schlichting  waves 
in  plane  and  curved  channel  flows.  Future  work  will 
study  the  nonlinear  stability  of  other  more  complex  2D 
and,  subsequently,  3D  channel  flows. 

Formulation  and  Numerical  Approach 

Governing  Equations 

The  temporal  evolution  of  3D  channel  flows  is 
governed  by  the  incompressible  Navier-Stokes 
equations.  These  equations  have  been  expressed  in  a 
general  nonorthogonal  coordinate  system  for  the 
purpose  of  this  study.  This  allows,  by  appropriate 
selection  of  a  coordinate  transformation,  the 
simulation  of  channel  flows  of  arbitrary  type  with  the 
restriction  that  both  the  geometry  and  the  flow 
solution  are  periodic  in  the  streamwise  and  spanwise 
directions.  It  is  also  assumed  for  present  purposes 
that  the  geometry  is  2D,  i.e.,  that  there  is  no  curvature 
in  the  spanwise  direction.  The  governing  equations  in 
general  coordinates  consist  of  the  contravariant 
components  of  the  incompressible  momentum 
equation,  and  the  continuity  equation.  These  are 


derived  by  Brandes-Duncan21  and  will  not  be  shown 
here.  The  dependent  variables  are  the  contravariant 
components  of  velocity,  along  with  the  pressure. 

Discretization  and  Solution  Method 

The  approach  to  the  spatial  discretization  of  the 
governing  equations  is  a  single-domain  spectral 
collocation  method,  and  uses  Fourier  discretization  in 
the  streamwise  and  spanwise  directions  and 
Chebyshev  discretization  in  the  normal  direction. 
Both  the  pressure  and  the  continuity  equation  are 
staggered  in  the  normal  direction,  and  Chebyshev 
interpolation  is  used  to  relate  values  at  the  staggered 
points  to  values  at  grid  points.  The  coordinate 
transformation  facilitates  the  Fourier  transform  of  the 
governing  equations  in  the  spanwise  direction,  and  this 
decouples  the  linear  parts  of  the  governing  equations 
from  each  other  in  the  spanwise  direction.  The 
equations  are  integrated  in  time  using  a  second-order 
accurate  time-centered  discretization.  A  fully-implicit 
iterative  method  is  used  to  solve  the  discrete  equations 
at  each  time  step,  to  the  desired  residual  convergence 
level.  All  of  the  nonlinear  terms  are  computed 
directly,  using  collocation  (matrix-multiplication) 
derivatives  and  are  lagged  in  the  iteration  procedure. 
For  details  of  the  formulation,  see  Brandes-Duncan21. 

High-Performance  Computing 

The  above  formulation  has  been  implemented  as 
part  of  an  object-oriented  framework  for  performing 
computational  fluid  dynamics  (CFD)  simulations, 
written  in  C++.  The  goal  of  the  object-oriented 
approach  is  to  achieve  maximum  code  re-use  between 
different  CFD  applications,  as  well  as  optimal 
computing  performance  across  various  scalar,  vector, 
and  parallel  computer  platforms.  Optimal  floating¬ 
point  performance  has  been  achieved  by  encapsulating 
machine-optimized  linear  algebra  routines  within  the 
implementation  of  CFD  objects,  while  retaining  a 
machine-independent  interface  for  using  those  objects. 
Details  of  the  object-oriented  approach  will  be 
presented  in  a  separate  paper.  The  direct  impact  of  the 
object-oriented  approach  on  the  current  study  is  that 
the  spectral  discretization  method  can  be  easily 
expressed  in  terms  of  lower-level  linear  algebra 
operations,  allowing  changes  to  the  discretization 
method,  or  to  the  linear  or  nonlinear  iteration  levels  of 
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the  solution  method,  to  be  made  easily  and 
independently  of  the  rest  of  the  code. 

Validation  of  Numerical  Approach 

Error  Analysis 

A  method  was  devised  for  computing  the  error  in 
the  discretization  scheme  by  prescribing  a  desired 
solution  and  analytically  computing  a  forcing  function, 
which  could  be  added  to  the  residual  such  that  the 
prescribed  solution  would  satisfy  the  modified 
governing  equations.  Prescribed  solutions  were 
selected  so  that  the  error  due  to  the  Fourier, 
Chebyshev,  and  temporal  discretizations  could  each  be 
isolated  and  examined  independendy.  For  example, 
the  prescribed  solution  for  examination  of  the  Fourier 
discretization  in  the  stream  wise  direction  is  shown  in 
Fig.  1.  Both  the  residual  level  and  the  solution  error 
after  10  time  steps  with  At=0.001  are  indications  of 
the  convergence  of  the  discretization  with  increasing 
number  of  points,  and  are  shown  in  Fig.  2.  Similar 
results  are  shown  for  the  Chebyshev  discretization  in 
Fig.  3.  For  both  the  Fourier  and  Chebyshev 
discretizations,  the  error  was  found  to  decrease  almost 
exponentially  with  the  number  of  points  and  machine- 
zero  convergence  was  achieved  with  less  than  52 
points.  The  error  analysis  was  performed  using  the 
residual  evaluation  and  implicit  iterative  solver  in  the 
CFD  code,  and  thus  saved  as  the  first  validation  of 
the  accuracy  and  correctness  of  the  code. 

Linear  Stability  Results 

The  second  validation  of  the  code  was  obtained  by 
computing  the  temporal  evolution  of  small  amplitude 
Tollmien-Schlichting  waves  in  2D  plane  and  curved 
channel  flows  and  by  comparing  the  resulting  growth 
rate  and  wavespeed  with  linear  stability  theory.  A 
spatial  resolution  of  8x48  (stream wise  x  normal)  grid 
cells  was  used  with  At=T/1000,  where  T  is  the  period 
of  oscillation  of  the  Tollmien-Schlichting  wave  from 
linear  stability  theory.  Shown  in  Table  1,  the  results 
of  this  study  demonstrate  comparison  with  linear 
stability  theory  with  agreement  up  to  4  or  5  decimal 
places. 


Nonlinear  Stability  Results 

Plane  Channel  Saturated  State 

The  first  stage  in  the  transition  to  turbulence  of 
plane  channel  flow  is  the  growth  of  the  primary 
Tollmien-Schlichting  wave.  Plane  channel  flow 
exhibits  subcritical  instability,  in  that  growth  of  finite- 
amplitude  instability  waves  occurs  for  Reynolds 
numbers  below  the  critical  Reynolds  number  predicted 
by  linear  stability  theory.  Early-time  results  for 
Reynolds  number  of  4000  and  streamwise 
wavenumber,  a,  of  1.25,  are  shown  in  Fig.  4.  The 
initial  conditions  are  Tollmien-Schlichting  waves  of 
different  amplitudes.  For  all  cases,  the  primary  mode 
initially  follows  a  linear  evolution,  while  the 
harmonics  of  the  primary  mode  experience  rapid 
growth  through  nonlinear  interaction  with  the  primary 
mode.  After  the  harmonic  modes  have  grown  to  a 
sufficiently  high  amplitude,  the  growth  of  the  primary 
mode  is  modified.  For  initial  amplitudes,  At ,  above  a 
critical  value  of  approximately  0.08,  the  primary  mode 
will  begin  to  grow,  and  will  continue  to  grow  over  a 
long  period  of  time,  until  it  reaches  a  saturated 
equilibrium  amplitude.  At  this  amplitude,  at 
approximately  At  =0.32,  the  Tollmien-Schilichting 
wave  will  persist  indefinitely  in  a  2D  simulation  as  a 
stable  flow  solution;  however,  this  wave  is  highly 
unstable  to  3D  disturbances,  and  will  break  down  into 
successively  complicated  states  in  a  physically 
realistic  3D  flow.  The  2D  state  achieved  after  20 
periods  of  the  Tollmien-Schlichting  wave  is  shown  in 
Figs.  5  and  6.  Though  the  wave  is  not  completely 
saturated  (saturation  would  require  hundreds  of 
Tollmien-Schlichting  periods),  the  .  growth  is 
sufficiently  small,  and  the  solution  is  representative  of 
the  final  state.  It  also  compares  well  with  the  time- 
asymptotic  solution  found  for  this  case  by  Orszag  and 
Patera3. 

Evolution  Equation  for  Plane  and  Curved 
Channels 

Theories  of  nonlinear  stability  examine  the 
nonlinear  interactions  that  occur  between  primary 
modes  active  in  the  flow  and  their  respective 
harmonics.  Weakly  nonlinear  stability  theory  uses  an 
asymptotic  expansion,  with  the  amplitude  of  the 
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primary  mode  as  a  small  parameter,  and  considers  the 
higher-order  effects  of  the  harmonics  on  the  evolution 
of  the  primary  mode;  for  example,  see  Singer  et  a/.13. 
In  these  approaches,  it  is  often  ascertained  that  the 
temporal  growth  of  the  primary  mode  is  governed  by 
an  ordinary  differential  equation  of  the  form: 


-Mk  = / '<*.> 


A,  dt 


(1) 


where  A,  is  the  magnitude  of  the  complex  amplitude  of 
the  primary  mode.  For  infinitesimal  amplitudes,  the 
evolution  equation  becomes  linear: 


—~Al  =  <o; 
Ax  dt  1  ; 


(2) 


Summary 

The  first  stage  of  development  and  validation  of  a 
numerical  tool  for  simulation  of  channel  flow  stability 
problems  was  completed.  Linear  stability  results, 
accurate  to  at  least  4  decimal  places,  were  obtained. 
Nonlinear  stability  results  were  compared  with  the 
Landau  equation  and  showed  that  weakly  nonlinear 
theory  applies  in  the  low  amplitude  range,  up  to  Aj 
=0(0.01),  and  the  evolution  equations  saturate  at 
higher  amplitudes.  Future  work  will  investigate  the 
nonlinear  stability  of  additional  types  of  channel 
flows,  and  will  extend  the  analysis  to  study  the 
instability  of  these  flows  to  3D  disturbances. 


with  growth  rate  <or .  For  small  but  finite  amplitudes, 
the  evolution  equation  is  the  Landau  equation22: 

-f  4*. m ‘ '  A‘  <3 

A  j  at 

where  the  Landau  coefficient,  /,  determines  the 
stability  of  the  solution  to  finite-amplitude 
disturbances.  In  particular,  if  u>,  and  l  are  negative, 
the  flow  exhibits  subcritical  instability  and  growth  will 
occur  for  values  of  A  ,  above  a  critical  value. 

The  evolution  equation,  Eq.  (1),  for  plane  and 
curved  channel  flows  was  investigated  numerically  by 
performing  multiple  simulations  with  varying 
amplitudes.  After  1  period  of  the  Tollmien- 
Schlichting  wave,  the  growth  rate  (1/A,)  dA,/dt  and  the 
square  of  the  amplitude  were  computed.  These  results 
were  compiled  for  various  cases  and  are  shown  in 
Figs.  7-9.  The  rate-axis  intercept  of  each  curve  is  the 
linear  growth  rate  ,  and  the  initial  slope  is  an 
estimate  of  the  negative  of  the  Landau  coefficient,  -/ . 
These  results  show  that  the  Landau  equation  is  valid 
only  for  small  amplitudes  up  to  approximaely  0.01. 
The  Landau  coefficients  for  five  different  cases  are 
shown  in  Table  2.  The  results  show  the  trajectory  of 
the  primary  mode  in  amplitude  squared  vs.  rate  space. 
This  trajectory  is  indeed  determined  by  an  ordinary 
differential  equation  of  the  form  of  Eq.  (1);  however, 
this  trajectory  is  not  polynomial  in  nature  and, 
therefore,  cannot  be  represented  with  a  finite  number 
of  polynomial  terms. 
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-16.573 

11.8 
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5772.22 
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1.02 

0.269215 
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0.269214 

0 

23.3 

Table  1.  Comparison  of  wavespeed,  co*  and  growth  rate,  ,  from  simulations  on  an  8x48  grid,  with  At=T/1000 

for  plane  and  curved  channel  flow.  Re  is  the  Reynolds  number,  Rc  is  the  radius  of  curvature  for  the 
curved  channel  cases  and  a  is  the  streamwise  wavenumber. 


Figure  1.  Prescribed  solution  for  testing  accuracy  of  the  Fourier  collocation  method  for 
the  plane  channel.  The  figure  shows  u,  v,  and  p  along  the  upper  wall  of  the 
channel. 
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Figure  2.  Residual  and  solution  error  for  the  Fourier  discretization  vs.  the  number  of 
points  used 


Figure  3.  Residual  and  solution  error  for  the  Chebyshev  discretization  vs.  the  number 
of  points  used 
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E,  /  E,(0) 


Plane  Channel,  Case  2,  Re=4,000,  a=1 .25 


Figure  4.  Growth  of  the  primary  disturbance  mode  for  one  period,  T,  for  plane  channel  flow  with  Re= 4,000, 
a=1.25  and  with  varying  initial  amplitude  At . 
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Figure  5.  Quasi-saturated  Tollmien-Schlichting  wave  in  a  plane  channel  with  Rs= 4,000,  <x=l  .25  after  20  periods 
of  the  wave:  contours  of  perturbations  in  (a)  u,  (b)  v,  (c)  p,  and  (d)  spanwise  vorticity,  as  well  as  total 
flow  quantities  (e)  u,  and  (f)  spanwise  vorticity. 


Figure  6.  Quasi-saturated  Tollmien-Schlichting  wave.  Perturbation  velocity  vectors 
and  streamlines. 
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Figure  7.  Curve  showing  the  numerically-computed  evolution  equation  for  plane 
channel  flow  at  the  critical  Reynolds  number,  Re=5772,  a=1.02. 


Figure  8.  Evoludon  equadon  for  plane  channel,  Re= 4,000,  a=1.25. 
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Figure  9.  Evolution  equation  for  curved  channel,  Re=4,00Q,  Rc= 1,000,  o=1250. 


:  Case 

A. 

1 

-1.40 

- 

2 

-8.56 

0.0077 

3 

-1.22 

- 

a 

-10.70 

0.0063 

5 

-1.16 

- 

Table  2.  Nonlinear  stability  results  for  the  plane  and  curved  channel  cases  shown  in 
Table  1.  The  Landau  coefficient,  /,  and  threshold  amplitude.  A, ,  were 
estimated  from  the  evolution  equation  results 
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Abstract 

The  purpose  of  this  paper  is  to  describe  the  discrete  problem  and  iterative  solution  method  for  spectral 
simulations  of  the  3-D  Navier-Stokes  equations  using  a  general  3-D  coordinate  mapping.  A  spectral  collocation 
discretization  is  proposed  along  with  a  multi-level  iterative  method,  and  has  two  unique  advantages.  First,  arbitrary 
coordinate  mappings,  which  introduce  non-constant  coefficients  into  the  governing  equations,  create  no  additional 
complication  in  the  spectral  collocation  approach,  but  are  not  feasible  for  spectral  Galerkin  and  Tau  methods. 
Second,  the  coupling  between  velocity  and  pressure  is  treated  implicitly,  allowing  either  fully-  or  semi-implicit 
schemes  to  be  used,  and  removing  the  need  for  the  artificial  velocity  boundary  conditions  required  with  fractional- 
step  methods.  The  performance  of  the  proposed  method  is  competitive  with  current  methods  in  both  memory 
utilization  and  operation  count.  Results  were  obtained  for  a  canonical  case,  curved  channel  flow,  to  test  the 
suggested  approach.  Accuracy  was  examined  via  the  predicted  growth  rate  of  a  Tollmien-Schlichting  wave.  Then, 
simulations  for  2-D  Dean  vortex  How  and  3-D  breakdown  of  a  saturated  Tollmien-Schlichting  were  obtained. 

Introduction 

Direct  numerical  simulation  (DNS)  of  problems  of  instability,  transition,  and  turbulence  in  flows  with 
simple  geometric  boundaries  such  as  channel,  boundary  layer,  and  pipe  flows  is  used  extensively  to  study  the 
complex  physical  phenomena  leading  to  turbulence.  Spectral  methods  are  typically  employed  for  solving  the  3-D 
incompressible  Navier-Stokes  equations,  resulting  in  an  accurate,  non-dissipative  discretization.  For  all  but  the 
simplest  geometry,  a  coordinate  mapping  must  be  used  for  the  computational  mesh.  The  resulting  form  of  the 
Navier-Stokes  equations,  however,  is  not  directly  amenable  to  spectral  Galerkin  or  Tau  methods,  due  to  the  presence 
of  non-constant  coefficients,  which  introduce  convolution  terms  into  the  wave-number  formulation  of  the  governing 
equations.  For  this  reason,  many  researchers  have  turned  to  spectral  element  methods  which  can  handle  more 
adequately  geometric  complexities — see,  for  example,  Guzman  and  Amon[l].  Several  studies  have  also  been 
reported  in  which  coordinate  mappings  have  been  used  along  with  single-  or  multi-domain  spectral  methods. 
Blodgett[2]  has  employed  a  multi-domain  spectral  collocation  method  in  a  2-D  simulation  of  leading-edge 
receptivity  using  a  streamfunction-vorticitv  formulation  and  a  Schwartz-Christoffel  (orthogonal)  transformation. 
Carlson  et  al.[ 3]  have  used  a  time-dependent  mapping  of  the  wall-normal  coordinate  in  3-D  DNS  in  a  channel  with 
an  emerging  obstacle  on  the  lower  wall.  They  used  a  Fourier  Galerkin  method  in  both  the  streamwise  and  spanwise 
directions,  and  a  Chebyshev  Tau  method  in  the  wall-normal  direction.  They  accomplished  this  by  treating  only  the 
Cartesian-like  constant  coefficient  terms  in  the  governing  equations  implicitly  and  updating  the  remaining  terms 
iteratively. 

The  present  study  proposes  an  extension  to  these  studies  by  allowing  for  an  arbitrary  3-D  coordinate 
mapping.  The  use  of  collocation,  rather  than  Galerkin  or  Tau,  derivatives  allows  the  implicit  part  of  the  algebraic 
system  to  include  coefficients  representing  non-constant  coordinate  metrics.  The  collocation  discretization  method  is 
described  below.  A  multi-level  iterative  approach  is  then  described  for  solving  the  resulting  nonlinear  algebraic 
system.  Finally,  2-D  and  3-D  simulation  results  are  presented. 
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The  Discrete  Problem 


Coordinate  Mapping 

The  3-D  incompressible  Navier-Stokes  equations  in  non-dimensional  form  are  solved  using  the  skew- 
symmetric  form  of  the  convective  terms[4],  and  are  written: 

^.  +  i(u.Vu  +  V  (ua))  =  -Vp  +  -J-V2u  ,  (I) 


V  •  u  =  0  .  (2) 

A  non-moving  general  coordinate  mapping  is  employed,  and  is  described  as: 

=Si(*|.*2>*3)»  for  >  =  (3) 

The  governing  equations,  Eqs.  (1)  and  (2),  are  written  in  generalized  coordinates  by  representing  vectors  and 
differential  operators  in  terms  of  covariant  and  contravariant  base  vectors  associated  with  the  above  transformation. 
See,  for  example,  Gal-Chen  and  SomervilIe[5],  for  a  more  thorough  discussion  of  generalized  coordinates.  The  base 
vectors  are  defined  as: 

dx : *  .  at. - 

Covariant:  e,  =  — ~i.-  ,  Contravariant:  e  =- — i.  ,  (4) 

dq,  dxj 

for  i=  1,2,3,  and  where  i;  are  the  Cartesian  base  vectors.  Automatic  summation  is  assumed  for  repeating  indices. 


The  covariant  and  contravariant  base  vectors  are  related  to  the  covariant  and  contravariant  forms  of  the  metric 
tensor,  respectively,  and  to  the  Jacobian  of  transformation  for  the  generalized  mapping  of  Eq.  (3),  which  will  be 
denoted  as  <Jg  .  These  relationships  are: 


Covariant:  g-  =  e,  •  e;  ,  Contravariant:  glJ  =  e'  *ey  , 


y[g  =  +[Det(g,y  )]"‘  =  +[Det(5'y  )]''/2  = 


Det 


3fe..§2.§3.y 


(5) 


Vectors  in  transformed  space  are  represented  alternatively  in  terms  of  Cartesian,  covariant,  or  contravariant 
components.  The  notation  used  is: 

Cartesian:  u  =  «,  if  ,  Covariant:  u  =  ufef  ,  Contravariant:  u  =  u'er.  .  (6) 

Derivatives  of  the  covariant  base  vectors  are  also  needed.  These  are  written  as: 


iL 

*7 


rme 

1  IJ  V/7 


(7) 


where  T^are  Christoffel  symbols  of  the  second  kind.  Contravariant  components  of  the  velocity  vector  are  used  in 

both  the  momentum  and  continuity  equation,  and  contravariant  components  of  the  momentum  equation  are  used  to 
represent  it  as  three  scalar  equations.  Appropriate  wall  boundary  conditions  on  velocity,  along  with  the  specification 
of  the  mean  pressure,  complete  the  continuous  problem.  See  Brandes-Duncan[6]  for  a  more  complete  derivation. 


Spatial  Discretization 

Following  Canute  et  ai [7],  the  continuous  system  is  spatially  discretized  using  a  spectral  collocation 
method  by  selecting  for  each  coordinate  direction  a  sequence  of  orthogonal  basis  functions.  For  a  staggered  grid 
arrangement,  the  velocity  is  represented  at  different  points  than  the  pressure,  and  sets  of  collocation  points  are 
selected  for  both  variables  in  each  coordinate  direction.  The  momentum  equation  is  solved  at  the  velocity  points, 
and  the  continuity  equation  is  solved  at  the  pressure  points.  Interpolation  on  the  spectral  basis  is  performed  to 
transfer  needed  quantities  between  the  two  grids.  A  concise  notation  for  the  spatially  discrete  system  is  to  denote  the 
spectral  collocation  derivative  along  the  c,  coordinate  as  the  matrix  operator  D\  interpolation  from  the  velocity  grid 
to  the  pressure  grid  as  the  matrix  operator  A*,  and  interpolation  from  the  pressure  grid  onto  the  velocity  grid  as  A0. 
Using  upper-case  letters  to  represent  the  discrete  values  of  the  velocity  and  pressure,  the  spatially  discrete 
momentum  equation  is  written  as: 


r  j  i  | 

•2— +yv'(U)  =  -C'(P)  +  —  L‘(U)  ,  where 
dt  Re 

JV'(U)  =  j (JgU  iUi  )+^U  JDjU ‘  +  T)kU  ‘ U *  . 

G'{P)  =  g^D^A" p)  ,  and 


(8) 
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The  spatially  discrete  continuity  equation  is: 

£>(U)  =  0  ,  where 


D(U)  =  A 


(9) 


Temporal  Discretization 

In  this  study,  both  fully-implicit  (FI)  and  semi-implicit  (SI)  second-order  time  integration  schemes  are 
considered.  For  both  schemes,  Crank-Nicolson  time-centered  integration  is  used  for  the  linear  terms.  For  the  FI 
scheme,  Crank-Nicolson  is  also  used  for  the  convective  term.  For  the  SI  method,  a  second-order  Adams-Bashforth 
method  is  used  for  the  convective  term,  which  results  in  a  linear  algebraic  system  at  each  time  step.  The  final 
discrete  form  of  the  momentum  equation  is  then  given  as: 


U'M+i  -U'J' 

2  Re L 

+y[c'(pn+1)+c'(P',)]=o  . 

where  for  the  FI  scheme: 

)  +  L'(U")] 

(10) 

NiM+\n  (U) = ijy  (u"+i ) + n‘  (ir )] 

and  for  the  SI  scheme: 

(ID 

Nij,+in(U)  =  ^[3N‘  (U")- N‘  (U""1)] 
The  continuity  equation  is  satisfied  by  the  velocity  at  every  time  step: 

(12) 

D(Un+l)  =  0  . 

(13) 

Initial  and  boundary  conditions  complete  the  specification  of  the  discrete  problem.  Since  velocity  points  are  located 
on  the  walls,  wall  boundary  conditions  on  velocity  are  implemented  by  replacing  the  discrete  equation  on  the  wall 
with  the  required  boundary  conditions.  Since  pressure  points  are  not  located  on  the  wail,  no  extraneous  wall 
boundary  conditions  on  pressure  are  needed.  To  specify  the  mean  value  of  pressure,  a  continuity  equation  is 
removed  at  an  arbitrary  location  and  replaced  with  the  condition  that  the  mean  value  of  pressure  be  zero. 


Iterative  Approach 


Outer  Iteration 

Low-memory  solution  strategies  for  solving  the  3-D  discrete  system,  Eqs.  (10)-(13),  require  the  use  ol 
multi-level  iterative  methods.  First,  the  nonlinear  algebraic  system  for  the  solution  at  the  n+1  time  level  in  the  FI 


method  is  solved  using  quasi-Newton’s  method;  where  the  Jacobian  matrix  is  obtained  by  linearizing  the  nonlinear 
terms  and  evaluating  them  at  the  prescribed  base  velocity.  The  quasi-Newton  iteration  level  is  written  as: 


'A  GlfSUl  \-R?\  Jjj/t+i.tfM-1  j  r5u| 

D  B  jSPJ  {-/?"'}  '  |pn+U,+l|  =  ■  (U) 

where  the  matrices  A,  D,  G,  and  B  represent  contributions  to  the  Jacobian  matrix  from  various  terms  in  Eqs.  (10)- 
(13).  Note  that  B  is  null  with  the  exception  of  a  single  row  which  enforces  the  mean-pressure  condition.  Also,  A 
contains  contributions  from  the  linearization  of  the  convective  term,  evaluated  at  some  base  velocity  solution.  The 
right-hand-side  vectors  Ri  and  Ri  represent  the  value  of  the  residual  of  the  momentum  equation  and  continuity 
equation,  respectively.  The  iterative  procedure  begins  (at  m= 0)  with  the  specification  of  an  initial  guess  for  the 
solution  at  the  n+ 1  time  level,  and  continues  until  some  norm  of  the  right-hand-side  is  below  the  convergence 
criterion.  For  the  SI  method,  the  above  method  is  also  used,  but  produces  a  converged  solution  in  one  iteration, 
since  the  algebraic  system  is  linear. 


Inner  Iteration 

An  inner  iteration  strategy  is  used  to  solve  the  linear  algebra  problem  resulting  from  the  quasi-Newton 
scheme.  With  the  linear  algebra  problem  abbreviated  as: 


LX  =  F  ,  where  X  = 


,  F  = 


(15) 


(16) 


a  straight-forward  iteration  technique  for  solving  this  system  is  written  as 

MSX  =  F-LXm  , 

X"'+l  =  Xm  +SX  , 

where  M  is  a  preconditioner  for,  or  rather,  the  some  easily  invertible  approximation  to,  the  coefficient  matrix  L.  The 
above  iterative  method  is  known  as  Richardson  iteration,  or  residual  correction.  Note  that  the  iterative  scheme 
requires  the  inversion  of  M,  rather  than  L,  and  that  L  is  only  needed  as  a  multiplicative  operator.  Other  more 
complicated  iterative  procedures  such  as  minimum  residual  methods  have  the  same  requirements  for  the  operators  M 
and  L.  For  the  inner  iteration  method,  the  preconditioner  M  is  chosen  as  an  approximate  lower-upper  factorization 
of  L,  of  the  form: 


M  = 


a  o  Ti 

G' 

’A 

AG' 

D  B-DgJLo 

I 

D 

B 

where 


07) 


M  differs  from  L  only  in  the  momentum  equation  pressure  coefficient.  The  following  procedure  is  used  to  solve  the 
linear  system  MY=C: 

AYf=C,  , 


(B-DG)>;  =C:-DYf  ,  (18) 

Y,  =  Y,*  -GK:  . 

The  approximate  factorization  Eq.  (17),  produces  a  solution  procedure,  Eq.  (18).  identical  to  that  used  in  the 
fractional  step  method,  as  shown  in  Perot[8]  and  Dukowicz  and  Dvinsky[9].  In  the  fractional  step  method,  however, 
M  replaces  the  original  linear  operator  L.  and  the  original  linear  system,  Eq.  (15),  is  not  satisfied.  Modified  velocity 
boundary  conditions  are  used  to  improve  the  accuracy  of  the  fractional  step  method  and  to  make  the  procedure  stable 
in  time[10].  These  can  be  implemented  directly  in  the  above  linear  system.  Thus,  the  outer-inner  iteration  method 
proposed  here  can  be  reduced  to  the  fractional  step  method  with  modifications  to  the  linear  system,  and  using  the  SI 
method  with  one  quasi-Newton  iteration. 


Velocity  and  Pressure  Solution  Methods 

In  order  to  solve  the  first  two  steps  of  Eq.  (18)  for  the  velocity  correction  Y  and  pressure  correction,  f:, 
respectively,  iterative  methods  are  again  utilized.  The  velocity  system  is  well-conditioned,  and  therefore  converges 
rapidly  with  any  reasonable  iterative  method.  A  preconditioned  residual  correction  method  as  in  Eq.  (16)  is  used, 
with  a  low-memory,  line-implicit,  preconditioner. 

The  pressure  system,  on  the  other  hand,  is  ill-conditioned,  and  constitutes  the  majority  of  the  computational 
effort  in  any  simulation.  A  preconditioned  GMRES  method  with  a  low-memory  line-  or  plane-implicit 
preconditioner  was  implemented,  and  was  effective  for  2-D  and  small  3-D  problems.  However,  the  number  of 
iterations  required  was  prohibitive  for  large  3-D  problems,  and  was  not  competitive  with  Galerkin  or  Tau  methods 


which  have  very  sparse  matrix  structures.  It  was  found  that  the  transformation  to  wave-number  space  could  be 
performed  on  the  linear  pressure  system  prior  to  solution  with  GMRES,  with  the  inverse  transformation  applied  to 
the  result.  The  transformed  pressure  system,  with  a  line-implicit  preconditioner,  converged  much  more  quickly  than 
the  original  one.  This  transformed  preconditioned  solver  for  pressure  allows  the  performance  of  the  present  scheme 
to  compare  reasonably  well  with  schemes  designed  for  simple  geometries  and  based  on  a  Galerkin  or  Tau  spectral 
formulation. 

Curved  Channel  Results 


Spectral  Discretization 


Solutions  using  the  3-D  spectral  collocation  formulation  with  the  outer-inner  iterative  method  are  shown  for 
a  periodic  curved  channel  problem.  Fourier  collocation  is  used  in  the  streamwise  (qO  and  span  wise  (^2)  directions. 
Chebyshev  collocation  is  used  in  the  wall-normal  (§3)  direction,  with  the  pressure  and  continuity  equation  staggered. 
The  computational  domain  is  then  defined  as: 


4,  e  [0.2tc]  ,  §2s[0,2ji]  .  e[-l.l]  . 


(19) 


The  metric  transformation  represents  the  curved  channel  geometry  using  scaled  polar  coordinates.  The  lengths  L\* 
L2,  and  L3  are  the  scaled  lengths  of  the  physical  domain,  where  L\  is  the  arc  length  of  the  channel  centerline.  The 
radius  of  curvature  of  the  channel  centerline  is  referred  to  as  R&  and  the  physical  domain  is  defined  as: 


0  € 


n 

2 


Ly  K 

2RC  2 


j±_ 

2  Rc 


y  e  [0,  ] 


r  e 


(20) 


The  typicallv-used  curvature  parameter,  r|,  is  related  to  Rc  throush: 


n=- 


Rr+  — 


.  Rr  = 


h.113. 

2  1-T) 


The  resulting  metrics  and  Christoffel  symbols  are  shown  in  Table  1. 
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Table  1.  Non-zero  metrics  for  the  curved  channel  mapping,  where  r  is  the  polar 


coordinate  radius. 


Accuracy  and  Stability  Issues 

The  nonlinear  algebraic  system.  Eqs.  (10)-(13),  was  solved  using  both  the  FI  and  SI  second-order  time 
integration  methods.  The  accuracy  of  these  two  methods  was  confirmed  by  computing  the  evolution  of  a  Tollmien- 
Schlichtirig  (TS)  wave  in  a  curved  channel,  with  Re =5 000,  and  r|=0.9802  (Rc=100),  the  streamwise  wavenumber  of 
the  TS  wave  was  a=!00.  The  amplitude  of  the  TS  wave  was  such  that  its  maximum  streamwise  velocity  was  \%  of 
the  mean  base  (low  streamwise  velocity.  At  this  amplitude  the  growth  rate  departs  slightly  from  that  predicted  by 
linear  stability  theory.  The  frequency  and  growth  rate  of  the  TS  wave  at  ;=0.1  are  shown  in  Table  2.  For  each  case, 
the  inner  iteration  was  fully  converged  (to  maximum  residual  <  10'12),  which  required  approximately  5  subiterations 
for  each  outer  iterations.  The  number  of  outer  iterations  required  for  the  FI  method  was  examined,  and  it  was 


observed  that  for  Ar=0.01.  3  outer  iterations  produced  a  fully-converged  residual.  Using  only  2  outer  iterations 
produced  a  method  which  was  less  accurate  in  time,  but  stable.  For  smaller  A t,  2  outer  iterations  were  sufficient  to 
converge  the  nonlinear  iteration.  The  FI  method  was  observed  to  be  unconditionally  stable,  while  the  SI  method 
required  a  time  step  on  the  order  of  Ar=0.001  to  insure  stability.  Other  validations  of  the  FI  method  are  described  in 
[6]  and  [11]. 


At 

Time-Marching 

Method 

COr 

0.01 

FI.  3  iterations 

0.375203 

3.29148  x  10'-' 

0.01 

FI.  2  iterations 

0.375210 

3.2908  i  x  10  ' 

0.001 

FI,  2  iterations 

0.375203 

3.29184  x  10-' 

0.001 

SI 

0.375203 

3.29185  x  10"' 

o.oooi 

FI,  2  iterations 

0.375203 

3.29189  x  10"' 

Table  2.  Frequency  ((Or),  and  growth  rate  (coO  values  at  /=0.1  for  a  curved  channel 
Tollmien-Schlichting  wave,  with  Re=5000,  r|=0.9802  (Rc=100),  a=100. 
The  inner  linear  iteration  was  fully-converged  for  all  cases.  FI  is  the 
fully-implicit  method,  SI  is  the  semi-implicit  method. 


Curved  Channel  Results 

Curved  channel  flow  has  been  examined  in  previous  studies  such  as  Finlay  et  ai [12],  to  demonstrate  the 
effect  of  streamwise  curvature  on  the  possible  routes  to  transition  to  turbulence.  In  this  study,  several  results  were 
obtained  for  curved  channel  (low,  and  are  listed  in  Table  3.  First,  2-D  results  in  the  cross-flow  plane  were  obtained 
for  saturated  Dean  vortex  solutions.  Curved  channel  flow  for  a  channel  with  at  least  moderate  curvature  exhibits  a 
streamwise  vortex  instability  called  the  Dean  vortex,  which  will  experience  linear  and  nonlinear  growth  above  the 
critical  Reynolds  number,  and  will  saturate  at  a  finite  amplitude.  For  Case  1  a  Reynolds  number  equal  to  the  critical 
Reynolds  number  of  1 14.26  was  used,  and  the  simulation  confirmed  that  the  Dean  vortex  mode  would  not  grow.  For 
Case  2,  with  a  higher  Reynolds  number  of  140.5,  the  Dean  vortex  grows  to  a  large-amplitude  steady-state  solution. 
Results  are  shown  in  Fig.  1.  For  a  higher  Reynolds  number  of  249.8  (Case  3),  a  second  pair  of  vortices  is  present  in 
the  steady-state  solution,  as  is  shown  in  Fig.  2. 


Case 

Re 

n,Rc 

a= 

2xc  Rc/Lt 

(3=2tc/L2 

Grid  nodes 
(N,xN2xN3) 

Maximum  t 
obtained 

1 

114.26 

0.975.  79 

0 

1.98 

1  x  9  x  49 

40 

2 

0.975,  79 

0 

2.5 

1  x  9  x  49 

500 

3 

0.975,  79 

0 

2.5 

1x9x49 

268 

4 

5000 

0.9802,  100 

100 

0 

9  x  1  x  49 

447 

5 

5000 

0.9802.  100 

100 

■m 

13  x  13  x  49 

527 

Table  3.  Curved  channel  flow  cases.  For  each  case  A/=0.01.  Values  of  0  for  a  and  (3  indicate  a  2- 


D  simulation  in  the  corresponding  plane. 

Secondly,  a  2-D  solution  of  the  saturated  TS  for  Case  4  was  obtained.  For  this  case  a  higher  Reynolds 
number  of  5000  was  used,  for  which  the  TS  wave  is  linearly  unstable,  but  grows  to  a  finite  amplitude  time- 
asymptotic  solution.  This  state  was  then  used  as  a  starting  solution  for  the  3-D  simulation  of  Case  5.  This  simulation 
was  performed  by  adding  a  streamwise  vortex  disturbance  to  the  saturated  TS  wave  solution.  Results  obtained  for  80 
characteristic  time  after  the  introduction  of  the  streamwise  vortex  mode  show  the  breakdown  of  the  TS  wave.  The 
energies  of  the  various  streamwise,  spanwise,  and  combination  modes  are  shown  in  Fig.  3,  with  the  typical  notation 
E(ki,k2)  for  each  energy  representing  the  krth  mode  in  the  streamwise  direction  and  the  k2-th  mode  in  the  spanwise 
direction.  The  increase  in  amplitude  of  even  the  highest  frequency  modes  is  typical  of  a  3-D  breakdown  to 
turbulence,  rather  than  bifurcation  to  a  3-D  equilibrium  state.  Further  3-D  results  were  not  available  at  the  time  of 
this  writing,  but  will  be  shown  in  the  conference  presentation. 


Conclusion 

The  present  study  described  a  new  approach  to  spectral  DNS  in  generalized  coordinates.  The  key  elements 
to  this  approach  are  the  collocation  discretization  and  the  fully-coupled  multi-level  iteration  scheme.  Future  work 
will  include  further  investigation  of  the  routes  to  turbulence  in  curved  channel  flows,  including  the  low  Reynolds 
number  phenomena  of  wavy  and  twisting  Dean  vortices,  and  the  higher  Reynolds  number  breakdown  of  TS  waves 
due  to  streamwise  vortex  disturbances.  Then,  similar  phenomena  will  be  investigated  in  the  presence  of  a  periodic 
array  of  surface  roughness  elements,  employing  a  3-D  non-orthogonal  grid  transformation  with  the  present  method. 


Figure  1.  Cross-section  of  the  velocity  field  for  the  Dean  vortex,  Case  2.  a)  Contours  of  the 
streamwise  velocity,  b)  Vectors  and  sample  streamlines  of  the  cross-section  velocity 
_ components.  The  solution  was  interpolated  onto  a  finer  grid  for  plotting  purposes. _ 


Figure  3.  Time  evolution  of  the  energies  in  various  modes  for  the  3-D  TS  wave,  streamwise  vortex 
interaction  simulation.  Case  5. 
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Abstract 

This  study  has  been  undertaken  to  analyze 
the  effect  of  compressibility  on  the  dynamic-stall 
phenomenon  by  accurately  simulating  the 
prevailing  mechanisms  of  the  formation  of  the  stall 
vortex.  Towards  this,  a  generalized  analysis  is 
developed  in  a  time-dependent  curvilinear 
coordinate  system.  It  is  implemented  with  flow- 
adaptive  gridding  and  arbitrarily  maneuvering  and 
deforming  bodies.  One  of  the  major  contributions 
of  this  analysis  is  the  elimination  of  the  boundary 
condition  of  zero  normal  pressure  gradient  at  solid 
surfaces,  which  breaks  down  in  regions  where 
separation  occurs  and  in  regions  of  high  surface 
curvature.  The  branch  cuts  are  treated  properly  by 
solving  the  complete  Navier-Stokes  equations  at 
these  interior  points.  The  preliminary  verification 
study  led  to  the  conclusion  that  computer 
resources  on  the  CRAY  Y-MP864  are  not 
adequate  for  the  analysis  being  undertaken.  This 
resulted  in  developing  an  object-oriented  linear 
algebra  class  library  for  high-performance 
computation  on  parallel  machines.  At  the  time  of 
the  conference,  the  verification  analysis  was  not 
completed,  and  as  such,  it  was  decided  to  provide 
some  information  towards  the  approach  used  in 
this  work. 

Introduction 

Dynamic  stall  occurs  when  an  airfoil  is 
pitched  rapidly  past  its  static  stall  angle.  The  lift 
as  well  as  the  drag  continues  to  increase,  to  a 


dramatic  extent,  past  the  maximum  possible  static 
value.  As  seen  in  the  experimental  results  of 
Jumper  et  al.  (1988)  in  Fig.  1,  as  the  angle  of 
attack  approaches  30°  during  the  pitch-up  cycle, 
there  is  a  dramatic  loss  in  lift,  and  a  corresponding 
and  sudden  ‘nose-down’  moment.  This  event, 
termed  dynamic  stall,  is  a  barrier  technology, 
limiting  the  operational  envelope  of 
supermaneuverable  aircraft  and  helicopters,  and 
the  performance  of  wind  turbines  and 
turbomachinery. 

As  is  understood  presently,  the  event  of 
dynamic  stall  for  incompressible  flow  is 
dominated  by  vortex  interactions  and  instabilities. 
The  instability  mechanisms  observed  in  two- 
dimensional  flows  will  be  described  below.  As  the 
airfoil  is  accelerated  from  rest  in  a  quiescent  flow, 
the  initial  instability  of  this  flow  occurs  when  the 
Reynolds  number  (Re)  exceeds  a  certain  relatively 
small  value  and  flow  separates  from  the  airfoil 
near  the  trailing  edge.  As  the  Reynolds  number  is 
increased  further,  this  initial  instability  can 
amplify  small  asymmetric  disturbances  in  the  flow 
and  an  asymmetric  instability  occurs  near  the 
trailing  edge,  eventually  developing  into  a  Karman 
vortex  street.  This  is  a  stable  periodic  flow  for 
low  enough  Re.  Due  to  the  periodic  nature  of  this 
flow,  the  specific  point  at  which  the  airfoil  is  given 
an  initial  angular  acceleration  to  begin  the  pitch-up 
is  extremely  important.  The  development  of  the 
resulting  flow  and  the  corresponding  lift,  drag  and 
moment  curves  are  observed  to  be  very  dependent 
on  when  the  pitch-up  is  begun. 
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For  flows  with  moderate  Re  values,  an 
instability  occurs  during  the  pitch-up  maneuver 
near  the  leading  edge  where  the  boundary-layer 
flow  separates,  eventually  intensifying  into  an 
attached  counter-rotating  upper-surface  vortex. 
This  vortex  is  the  cause  of  the  increased 
performance  of  the  airfoil,  in  terms  of  its  lift. 
Momentarily,  it  remains  attached  over  the  center 
of  lift,  resulting  in  very  little  change  in  moment 
applied  to  the  airfoil.  An  underlying  paired  vortex 
forms  on  the  airfoil  surface.  As  the  maneuver 
continues,  an  additional  pair  of  vortices  forms 
below  this  one.  When  the  second  vortex  of  this 
new  pair  is  formed,  the  initial  secondary  vortex 
(paired  with  the  dynamic-stall  vortex)  erupts  out 
into  the  nearly-inviscid  vortex  flow,  cutting  off  the 
feeding  shear  layer  emanating  from  the  leading 
edge,  and  the  dynamic  stall  vortex  evolves  as 
shown  in  Fig.  2.  The  shear  layer  undergoes  a 
Kelvin-Helmholtz  instability,  rolling  up  into  the 
detached  vortex  which  convects  downstream  over 
the  airfoil,  resulting  in  the  large  ‘nose-down’ 
moment  on  the  airfoil.  The  dynamic-stall 
phenomenon  in  incompressible  flow  has  been 
simulated  very  satisfactorily  by  K.  Ghia,  Yang, 
Osswald  and  U.  Ghia  (1992)  and  a  suitable 
modulated  suction/injection  control  law  was 
devised  by  Yang,  K.  Ghia,  U.  Ghia  and  Osswald 
(1993). 

The  presence  of  compressibility,  three- 
dimensionality,  transition  and  turbulence,  as  well 
as  a  number  of  other  factors,  can  significantly 
affect  the  resulting  flow  evolution. 
Compressibility  can  have  significant  effects  on  the 
sequence  of  instabilities  described  above.  For 
flow  with  a  relatively  low  free-stream  Mach 
number,  interaction  with  a  rapidly  maneuvering 
airfoil  locally  causes  the  flow  to  accelerate  to  a 
Mach  number  of  two  and  higher,  resulting  in 
significant  compressibility  effects,  specifically  the 
production  of  additional  vorticity.  Maneuvers 
performed  in  supersonic  free-stream  flows  have 
even  more  complex  flow  physics.  In  either  case, 
as  shown  in  Fig.  3,  shocklets  formed  from  local 
accelerations  of  the  flow  around  vortices  and  wavy 
unstable  shear  layers  can  significantly  affect  the 
resulting  flow  structure  according  to  the  results  of 
Chandreshekar  et  al.  (1995).  These  shocklets  exist 
over  only  a  small  range  of  angles  of  attack  during 


the  maneuver.  Small,  temporally  evolving  shocks 
forming  in  the  presence  of  vortex  interactions  and 
instabilities  present  challenging  flow  physics  to  be 
simulated  by  the  computational  fluid  dynamicist. 

In  order  to  accurately  predict  such  a  flow 
numerically,  the  instability  mechanisms  present  in 
the  incompressible  flow,  and  any  additional 
mechanisms  present  in  the  compressible  flow, 
must  be  accurately  resolved  both  spatially  and 
temporally.  Compressible  phenomena,  such  as 
shocklets,  must  be  accurately  predicted  in  order  to 
understand  their  interaction  with  vortices  and  their 
effect  on  instability  mechanisms.  The  effects  of 
assumptions  made  in  an  analysis  on  the  accuracy 
of  resolving  the  instability  mechanisms  must  be 
addressed.  Because  of  the  sensitive  nature  of 
instabilities  to  disturbances,  errors  made  in  the 
analysis  may  be  amplified  to  a  point  where  they 
affect  the  resulting  flow  physics.  The  standard 
zero  normal  pressure  gradient  boundary  condition 
used  on  solid  surfaces  in  compressible  viscous 
flow  calculations  such  as  those  performed  by 
Visbal  (1990)  and  Choudhuri  and  Knight  (1995) 
may  have  a  significant  effect  on  the  resulting  flow 
evolution,  because  the  eruption  event  described 
above  is  characterized  by  a  region  of  high  normal 
pressure  gradient  relative  to  the  wall.  Since  the 
unsteadiness  is  of  prime  interest  in  dynamic  stall 
flows,  time-accuracy  of  the  analysis  is  of 
fundamental  importance. 

The  objective  of  this  study  is  to:  i)  develop 
an  analysis  that  is  applicable  to  dynamic  motion 
CFD  problems,  ii)  to  use  object-oriented 
programming  techniques  in  the  development  of  the 
analysis  to  facilitate  the  use  of  complex  CFD 
techniques  (such  as  flow-adaptive  grids,  Chimera 
grids,  etc.)  on  ever  changing  high-performance 
massively  parallel  computer  platforms,  and  iii)  to 
develop  simulation  results  for  dynamic-stall 
phenomena  under  a  variety  of  flow  conditions  to 
help  in  the  understanding  of  the  role 
compressibility  plays  in  the  phenomenon  of 
dynamic  stall. 

Underlying  the  entire  analysis  are  the 
physical  issues  of  how  to  accurately  capture 
instabilities  and  avoid  error  growth,  as  well  as  the 
computational  issues  of  optimization  and 


efficiency.  The  computational  issues  should  not 
be  allowed  to  affect  the  analysis  in  such  a  way  as 
to  .  lower  confidence  in  the  resulting  solutions 
through  the  addition  of  errors. 


Analysis  and  Discussion 

The  analysis  begins  with  the  first 
principles  of  classical  field  theory  applied  to  a 
continuum  fluid.  The  integral  equations  are 
transformed  into  a  set  of  non-linear  partial 
differential  equations  written  in  terms  of  an 
arbitrarily  moving  coordinate  system. 
Unsteadiness  of  the  coordinate  system  as  well  as 
that  of  the  bodies  in  the  flow  are  accounted  for  in 
this  manner.  Cross  derivatives  are  retained  in  the 
resulting  differential  equations  to  ensure  accuracy 
of  the  solutions  of  these  equations.  The  equations 
are  neither  modified  by  explicitly  adding  damping 
terms  nor  by  using  a  turbulence  model.  Both  of 
these  are  inconsistent  with  the  very  equations  we 
are  trying  to  solve.  A  consistent  set  of  boundary 
conditions  have  been  developed  which  removes 
the  requirement  for  a  pressure  boundary  condition 
on  solid  surfaces.  Consistent  far-field  conditions 
are  used,  unlike  the  zerolh-order  extrapolation  or 
characteristics-based  boundary  conditions  com¬ 
monly  used  . 

The  differential  equations  are  discretized, 
resulting  in  a  non-linear  algebraic  equation  set 
whose  solution  approaches  the  continuum  solution 
at  a  known  rate.  Based  on  the  results  of  a  model 
problem,  it  can  be  said  that,  as  the  grid  is  refined, 
the  physics  of  the  discrete  equation  set  does 
approach  that  of  the  continuum  equations.  The 
discrete  equations  are  solved  at  all  points  within 
the  flow ,  including  pseudo-boundaries  introduced 
by  computational  topologies  such  as  branch  cuts 
and  multi-block/Chimera-grid  boundaries.  Time 
accuracy  is  assured  by  complete  convergence  of 
the  residual  of  the  non-linear  equation  set, 
guaranteeing  its  solution,  and  through  the  use  of  a 
fully  implicit  time  discretization.  Errors 
introduced  through  incomplete  convergence  of  the 
equations  may  adversely  affect  the  flow  physics  of 
the  discrete  equation  set  and,  thus,  the  instability 
events  that  are  being  captured. 


Numerical  Jacobians  are  used  to 
efficiently  generate  the  resulting  linearized 
coefficient  matrix.  A  fully  direct  method  is  not 
feasible  for  solving  the  resulting  very  large  system 
of  equations,  even  when  the  sparse  nature  of  the 
system  is  considered.  Thus,  an  iterative  method  is 
used  to  drive  the  residual  of  the  equation  set  to 
zero  at  each  time  step. 

Computational  optimization  pervades 
every  aspect  of  a  simulation  of  a  challenging  flow 
problem.  To  lower  memory  requirements,  the 
number  of  operations  and  numbers  of  iterations 
performed,  flow-adaptive  grid  techniques  are  used 
to  achieve  near-optimal  grid-point  placement  and, 
more  importantly,  to  accurately  capture  unsteady 
flow  phenomena.  Multi-grid  techniques  have  been 
programmed  and  will  be  used,  if  necessary,  to 
accelerate  the  convergence  of  the  unsteady 
equations  at  each  time  step.  If  the  problem 
requires  the  use  of  a  direct  method  for 
convergence-related  issues,  then  a  multi-block 
‘divide  and  conquer’  approach  is  planned  to  break 
the  problem  into  smaller  pieces.  In  addition, 
underlying  all  numerical  computations  is  the  issue 
of  matching  algorithms  to  the  machine  architecture 
in  use.  The  use  of  massively  parallel  machines  is  a 
necessity  for  computing  dynamic  motion  CFD 
problems.  Because  of  the  need  to  change 
algorithms  when  switching  architectures,  code 
design  becomes  important,  to  reduce  the  time 
spent  porting  between  various  machines. 

Object-oriented  programming  (OOP) 
techniques  have,  therefore,  been  used  to  address 
many  of  these  optimization  issues.  Code 
complexity  increases  dramatically  as  complex 
CFD  techniques,  such  as  flow-adaptive  and  multi¬ 
block  grids,  are  added  to  an  analysis.  Advanced 
programming  techniques  are  necessary  to  manage 
such  complexity  in  an  organized  manner.  The 
interaction  between  various  codes,  such  as  a  grid 
generator  and  a  flow  solver,  must  be  expressed  in 
the  program  itself,  in  order  to  allow  effective 
communication.  With  the  use  of  OOP  languages, 
portability  issues  are  more  directly  addressed  and 
more  adequately  resolved.  The  use  of  an  object- 
oriented  language  actually  facilitates  efficient 
programming  rather  than  hindering  it,  since  the 
way  an  object  is  used  is  decoupled  from  the 


particular  way  it  was  implemented  on  a  given 
machine.  Specific  optimized  algorithms  can  be 
used  where  appropriate,  without  complicating 
other  parts  of  a  code.  A  good  example  of  this  is 
linear  algebra  programming,  which  forms  the  basis 
of  all  CFD  codes.  A  particular  efficient 
implementation  of  a  matrix  operation  can  be  used, 
where  one  was  not  before,  without  changing  the 
way  the  CFD  program  loads  up  its  coefficient 
matrix.  An  example  of  the  high  performance 
attainable  on  a  CRAY  Y-MP  with  C++  is  given  in 
Fig.  4.  Because  of  the  design  of  the  linear  algebra 
objects  used  by  the  CFD  code,  this  optimized 
block  matrix  routine  was  used  without  changing  a 
single  line  of  code  in  the  CFD  program. 

The  use  of  user-defined  types  allows  the 
program  to  be  written  in  terms  of  objects  closer  to 
the  problem  at  hand  and  allows  the  compiler  to 
know  when,  for  example,  a  grid  was  misused  and 
mark  the  line  with  an  error  flag.  In  a  procedural 
language,  such  as  FORTRAN,  the  compiler  has  no 
way  of  knowing  if  a  programmer  really  meant  to 
pass  down  to  a  routine  the  x-coordinates  of  a  grid 
first,  then  the  y-coordinates,  or  vice-versa.  All  it 
knows  is  that  it  was  expecting  two  arrays,  and  is 
ignorant  of  what  they  were  supposed  to  be 
representing. 


Summary 

The  compressible  flow  past  a  maneuvering 
airfoil  is  being  studied.  Currently,  issues  related  to 
boundary  conditions,  including  proper  treatment  of 
branch  cuts  and  multi-block/Chimera  boundaries 
are  being  addressed.  Numerical  analysis  issues  for 
solving  large-scale  iterative  problems  are  being 
examined.  Finally,  object-oriented  programming 


is  being  used  to  facilitate  the  incorporation  of 
complicated  CFD  techniques  into  flow  codes  as 
well  as  to  address  issues  of  software  portability  on 
high-performance  supercomputers.  Results  of  this 
effort,  which  were  not  ready  at  the  time  of  the 
Conference,  will  be  disseminated  shortly  in  the 
form  of  a  report. 
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a)  Initial  leading  edge  separation 


b)  Onset  of  interaction,  sharp  local  pressure  change; 


c)  Strong  interaction  d)  inviscid  Interaction 


Figure  2. -Vortex  Structure  in  Sequence  of  Events  for  Dynamic  Stall  (K.  Ghia,  et  al„  1992) 


Figure  3. -Compressibility  Effects  on  Dynamic  Stall  t  Chandresekhara.  et  al..  1993) 


Figure  4. -Cray  Y-MP  Performance  for  C++  Block  Matrix  Inverse 


